2014-02-16 21:01:49 8 Comments

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.

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:

### Related Questions

#### Sponsored Content

#### 48 Answered Questions

#### 1 Answered Questions

### [SOLVED] Finding holes in 2d point cloud

**2018-05-17 12:46:22****bestyasser****387**View**0**Score**1**Answer- Tags: python geometry 2d computational-geometry point-clouds

#### 4 Answered Questions

### [SOLVED] Find point which sum of distances to set of other points is minimal

**2011-01-05 22:27:04****Pawel Markowski****8212**View**10**Score**4**Answer- Tags: algorithm optimization geometry gis

#### 3 Answered Questions

### [SOLVED] Given a vector of points (possibly out of order), find polygon (not convex hull)

**2011-09-13 21:04:13****Everaldo Aguiar****5687**View**8**Score**3**Answer- Tags: c++ algorithm geometry computational-geometry

#### 3 Answered Questions

### [SOLVED] Intersection of a line and a concave polygon 3D

**2014-06-13 14:57:36****DOFHandler****1105**View**3**Score**3**Answer- Tags: performance algorithm computational-geometry

#### 1 Answered Questions

### Set operations (union and intersection) on simple polygons in 2D

**2014-04-18 08:30:50****Peter K.****1215**View**1**Score**1**Answer- Tags: algorithm computational-geometry

#### 1 Answered Questions

### [SOLVED] Union of two polygons defined by paths

**2012-06-19 14:26:01****Charles****1361**View**0**Score**1**Answer- Tags: algorithm language-agnostic geometry computational-geometry

#### 4 Answered Questions

### [SOLVED] Polygon enclosing a set of points

**2009-05-06 09:57:53****Laurent K****31746**View**31**Score**4**Answer- Tags: algorithm language-agnostic geometry polygon

## 8 comments

## @Luke Hutchison 2020-10-11 00:45:46

`r`

to be the smallest radius that should lead to the detection of a hole (or half the smallest diameter of a hole in some direction).`x`

and`y`

in a doubly-nested loop, Perform a grid search over your image. The step size`s`

should be half of`r`

, or even much smaller, depending on the tradeoff you want to make in the ability to detect small holes (and to detect all points on the boundary of a hole) vs. the processing time.`(x, y)`

, find the nearest neighbor point`(nx, ny)`

using the k-d tree.`r`

, then the grid point is inside a hole, and the nearest neighbor is on the boundary of the hole. Mark these points as such.You can then use the Union-Find datastructure (or just use a depth-first flood fill algorithm) to find all the connected component grid points that are inside holes. You can then find the outline of the hole by tracing around the outline of the grid points that comprise the holes, and stringing the nearest neighbor points together in order (you might have to tweak this order a bit to produce a non-self-intersecting polygon).

The nearest neighbor algorithm might miss some of the points that are arguably on the boundary of the hole, but this will be less likely the smaller the value of

`s`

.## @Spektre 2014-02-19 14:51:44

what about some

bitmap+vectorapproach like this:obtain bounding box of point cloud area coverageDo this if it is not already known. It should be simple

`O(N)`

cycle through all points.create`map[N][N]`

of the areaIt 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 sizeN is also the accuracyof the output andmust be bigger then average point distance !!!so eachcellinside`map[][]`

covers area with at least one point (if not in hole area).compute data density for each cell of`map[][]`

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)`

create list of unused area`(map[][].cnt==0)`

or`(map[][].cnt<=treshold)`

I do it by Horizontal and Vertical lines for simplicity

segmentate outputJust group lines of the same hole together (intersecting ones ... vector approach) and also can be done in bullet

#4by flood fill (bitmap approach)polygonize outputTake all edge points of

H,V linesof 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.My source code for this approach:Just ignore my

`glview2D`

stuff (it is my gfx render engine for geometry)`view.pnt[]`

is dynamic list of your points (generated by random)`view.lin[]`

is dynamic list outputH,V linesfor visualization only`lin[]`

is your lines outputThis is 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

but for most of the code it is

`O(N)`

and for hole search/segmentation`~O((M^2)+(U^2))`

where:`N`

is point count`M`

is map grid size`U`

isH,V linescount dependent on holes ...`M << N, U << M*M`

as you can see for

`3783`

points`30x30`

grid on the image above it took almost`9ms`

on my setup[Edit1] played with vector polygonize a littlefor simple holes is fine but for more complicated ones there are some hick-ups yet

[Edit2] finally got a little time for this so here is it:This is simple class for hole/polygon search in more pleasant/manageable form:

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<int> x;`

is the same as`int x[];`

`x.add();`

add empty item to x`x.add(a);`

add a item to x`x.reset()`

clears the array`x.allocate(size)`

preallocate space to avoid reallocations on the run which is slow`x.num`

is number of items in x[] ... used size in itemsin the original code are only static arrays so if you are confused check with it instead.

Now how to use it:where

`view.pnt[]`

is list of input points and inside it:`view.pnt[i].p0.p[ 2 ]= { x,y }`

Output is in

`h.lin[]`

and`lin_i0`

where:`h.lin[i] i= < 0,lin_i0 )`

are the inside H,V lines`h.lin[i] i= < lin_i0,h.lin.num )`

are the perimeterThe 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:`view.lin[]`

has members:`p0,p1,`

which are points as`view.pnt[]`

and`col`

which is colorI saw only one issue with this when holes are too small

`(diameter < 3 cells)`

otherwise is OK[edit4] reordering perimeter linesto do that just instead of this:

do this:

[Edit5] How polygonize inside`holes::holes_end`

worksAs input for this you need the list of all

H,V lines`lin[]`

segmentated/grouped/sorted by hole and the density map`map[][]`

.loop through all holesloop through all H,V lines of processed holeCreate 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.delete all non border points from listSo remove all points that have no contact with non-hole area by looking into 4 neighbors in the density

`map[][]`

do connected components analysis on the pointsset`used=0; p0=-1; p1=-1;`

for all points in`pnt[]`

listconnect points with`distance=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 soaddthe connection info to booth of them (use`p0`

or`p1`

index which ever is unused`(-1)`

) and increment usage of both points.try to connect points with`distance=sqrt(2)`

is almost the same as

#2except 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.try to connect closest pointsagain is almost the same as

#2,#3but select the closest point instead and also avoid closed loops.form polygon from`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[]`

.## @Spektre 2015-12-02 14:46:21

@sendreams what do you mean by that? You want the hole perimeter polygon (that is what I do in this answer) or something else (like fill the non hole area... you need triangulation for that as the result is polygon with holes)? My reference to bounding box is rectangular area covering all points from dataset

## @Spektre 2015-12-04 08:30:24

@sendreams added [edit5] for you it is mainly compilation of the comments from the code itself + some additional explaining.

## @sendreams 2015-12-04 08:52:30

thanks for your explain. can you tell me some link about open source or libs?

## @Spektre 2015-12-04 08:59:46

@sendreams no because I code almost everything myself because lot of my work is really crazy and alien to common math/programming stuff so most things I do are not yet developed/derived/coded or if they are it is not usable on standard HW or time I would spend learning to use or change it to suit my needs would be far bigger then code it on my own... not to mention legal problems ... And the rest like this is so easy that I do it for fun and to ease my mind from the hard stuff anyway.

## @AndyTheCornbread 2018-07-12 18:04:55

I went ahead and accepted this as the solution. I ended up getting really sick shortly after I posted this and am finally recovering a bit now four years later. I haven't tested this solution because I stopped working on this while I was sick but it looks like it would do what I wanted and or get me a long ways down the path to accomplishing what I wanted. Sorry it took me so long to get logged back in here and look at this.

## @Nuclearman 2014-02-18 10:05:05

You'd probably be better off using the Delaunay triangulation to find the Gabriel graph. You then angle sort the Gabriel graph and do circular walks to generate a list of convex polygons. You can then sort those polygons by area. You'd be interested in the ones with the largest area.

It'll also be more efficient to modify the angle sorted graph that that you can follow the path from A to B, and see what is next either clockwise or counter clockwise (from the angle sort). A dictonary of dictionaries might be useful, that is defined something like "graph[A][B] = (clockwise, counterclockwise)". An example algorithm (python) using the dictionary of dictionaries approach.

It may also be useful to combine with Ante's suggestion as well.

## @Yves Daoust 2014-02-17 13:35:54

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 restore 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

## @piXelicidio 2014-02-16 22:10:45

This is my enthusiast non-scientific solution:

1 - Scan all the 2D area with minimum predefined step (dx, dy). For each step coord find the bigger circle that could fit without any point inside. Discard all the circles with radius less than a predefined size.

2 - Now find all groups of colliding circles, easy test of distance and radius, store and group in separated lists. (Ask, if you want more details about how to grouping them, is really easy )

3 - Find the concave bounding polygon for each group of circles, very similar to the algorithm to find the convex polygon around a group of points you already wrote, and your last question angles between vectors was related.

Notes

Optimization tips: Before step 1, you can store all points in a grid|matrix so distance calculation are simplified and limited to near grid squares of the given circle radius.

Precision: You gain more precision for smaller values of scan step and minimum circle radius allowed.

Not tested by myself but, I'm sure it works. Good luck!

## @Ante 2014-02-16 22:08:15

Delauney triangulation can help. It has property that no input point is inside the circumcircle of any triangle in triangulation. Because of that, hole boundary points will be connected by larger/wider triangles covering that hole. In your cases triangulation will have lot of triangles of similar size, and some triangles of larger size that cover holes. Probably it is enough to filter larger ones and connect them to find a hole.

## @AndyTheCornbread 2014-02-16 23:59:33

This is a solid idea. I have done Delauney Triangulation and Voronoi diagrams before for terrain modeling. I think using this and amit's idea of deviations from the mean as the parameter for adjustment I might be able to come up with something workable. It might be a few days before I get a chance to tackle this but I'll post back with how well this idea works.

## @Ante 2014-02-17 08:03:15

Consequence of this approach is also finding of exterior boundary. It is always a question when to stop finding holes, or what is the largest hole to consider. I suppose that radius of circumcircle or triangle longest edge is good measure to check is triangle a hole. Probably most of values will be similar, than mean will give value that is not a hole. If it is a case, than larger value (few times) can be used as a threshold.

## @sendreams 2015-12-02 13:22:40

it can be used to find concave hull, i think it's no good to find hole.

## @amit 2014-02-16 21:31:30

Here is a thought:

`x`

find the distance`d(x,y)`

(where`y`

is the closest neighbor to`x`

). Define`f(x)=d(x,y)`

as above.`f(x)`

.`f`

values are very far from the mean, far by at least \alpha standard deviations. (\alpha is a parameter for the algorithm).## @Spektre 2014-02-20 07:54:02

if the points are not spatially sorted and are a lot of them then this can be a very very slow approach (just first bullet is O(N^2) in that case) or am I wrong and missing something?

## @amit 2014-02-20 08:22:16

@Spektre Closest neighbor can be done in

`O(nlogn)`

using k-d trees. All the rest bullets are`O(n)`

## @Kyle 2014-02-16 21:31:06

I don't know of any algorithm off the top of my head, but one thing you might try (and this is my first impulse here) is something similar to how density is computed in meshfree methods like smoothed-particle hydrodynamics.

If you could compute a density value for any point in space (and not just at the sample points you're given) you could define the boundaries of the holes as level-set curves of the density function. I.e. the holes are where the density drops below some threshold (that you'd likely allow the user to configure). You could find the boundaries using something like marching squares.

If you'd like any clarification on how these sorts of density interpolation functions work, I can provide that (to the best of my ability and knowledge) if you like.

I don't know how well this would actually work, but hopefully it'll give you some direction.