I have 2D-points, which boundaries depict a polygon. The points are at integer positions (no fractions possible).

Now I want to build an enumerator in C# that returns all points that are inside the polygon. Implementing the enumerator is not the issue, but how do I efficiently get all the points? The polygon can be quite big.

  • Is polygon always convex or can be non-convex? – Kostya Regent Jan 5 '15 at 12:49
  • Good question, frankly do not know, yet....the polygons represent areas, which are randomly being generated. If convex would benefit in a simpler algorithm, I can live with it, I guess – Marco Jan 5 '15 at 13:00
  • There are at least 2 different solutions answered, but I want to add one more comment. If you would have problem like "how many points with integer coordinates is inside polygon" then you can solve it by dividing polygon to square triangles and rectangles. In this case you can easy calculate count of points (using area formulas) without iterating. This algorithm depends on count of polygon vertices, not on count of points inside polygon. It works with convex and non-convex polygons, but in non-convex case count of points in some triangle should be negative. – Kostya Regent Jan 5 '15 at 13:52
  • Ok, but I need the coordinates of the points. So, actually the polygon is in OO-terms a container and the 'points' inside are cells. Cells are childs of the container (in this case a polygon). But with the answers I can continue..thx – Marco Jan 5 '15 at 13:55
up vote 7 down vote accepted

What you are seeking is a solution to the "point-in-polygon" problem. It is described in Wikipedia here:


          Wikipedia image


You can find C code following the links to Computational Geometry in C, or at many other locations found by searching for "point-in-polygon" code.

  • Joseph, I do not search if a point is in a polygon, I want all points, which are inside a polygon. My question is what an effective algorithm would look like to tackle this problem – Marco Jan 5 '15 at 13:30
  • 2
    For each point, check if it is in the polygon. There is no simpler method, unless your points are organized in some nonrandom fashion that you did not specify. – Joseph O'Rourke Jan 5 '15 at 13:38
  • 2
    @Hawk66, you can improve this solution for your purposes. Make test for each Y coordinate polygon has and iterate by points when count of intersected lines is odd (odd means inside). – Kostya Regent Jan 5 '15 at 13:41
  • Ok, thx for the clarification – Marco Jan 5 '15 at 13:44
  • 2
    There are many point-in-poly algorithms listed here in C. tog.acm.org/resources/GraphicsGems/gemsiv/ptpoly_haines/… – Syntac_ Jan 5 '15 at 13:56

You will need to process all the points at least once so if this check is done only once there isn't much you can do to speed up the test other than brute-forcing it using parallelism.

If the test is going to be run multiple times there are ways to pre-calculate tables to help, such as a grid of cells marked as [definitely inside (green in image), outside (red), and maybe (blue)], only points that fall in the blue cells need to be tested against the polygon.

enter image description here

But if it's run only once with just few points relative to the grid or too many "maybe" cells then building this table is likely to take more time than just brute-force testing your points.

You need to have many more points than cells (low resolution grid) for this to turn a "profit", and most cells must end up as definite inside or outside (high resolution grid), so finding a balance that speed things up between those two factors can be difficult.

It's easier and more consistent in terms of gains to just split the number of points between multiple threads and do a simple odd-even intersection test. GPGPU processing such as OpenCL can give you an easy and huge speed up for this kind of test simply by using a massive number of cores.

  • If the points represent the corners of the grid's squares, it looks like there might be points that are inside your star but not depicted in green. They would be the corners of squares, so not all four points per square would be inside the star. The end solution might require half-pixel offsets or some such, such that any square whose center was inside the polygon was colored a certain way. – Seth Battin Jan 5 '15 at 19:02
  • 1
    @SethBattin It is my understanding that the blocky shapes in the diagram are NOT the unit-squares themselves, but larger regions. All unit squares (or even floating point coordinates) within each green block are inside the polygon, all unit squares within a red block are outside, and only unit squares within a blue block require a full test. The challenge in this approach is finding the right balance between evaluation speed (at small region size) and memory efficiency and initialisation speed (at large region size) for the particular use case. – DeveloperInDevelopment Jan 5 '15 at 20:58

Every polygon can be represented as lines[L] with interception points[P].For each line you calculate weather it is above or below a point Pn(n is the index of array P) which is not part of the line Lm.When you measure all the lines you save the result in one dimensional bool array V[number_lines].When you have a random point A.Then you make for loop for each line and you make the same check if point A is below or above line Li.As you do it you compare the result to V[i] and if its not the same you return true if point A is over the line Li else you return false.Once the loop though all lines from [L] you can return true.As you do this you can know if a point is inside or outside an polygon easly which means that you can loop though all A points in the polygon bound box Box[Px(min),Py(min),Px(max),Py(max)].To detect weather point A is above or below line L you must turn line L into function f(x)=y=ax+b by using the two points representing the line as x=P1.x y=P1.y and x=P2.x x=P2.x including them into system.Once you have calculated the function f(x) you must simply test if f(A.x) ><= A.y . If f(A.x)>A.y then the line L is above point A.If =0 then its on the line else its below it. It might look a bit complicated, but the most of those calculations are instant and can be calculated only once. As you store V,L,P arrays in the polygon represented class.

PS: If line L with points P1 and P2 is such so P1.x=P2.x then you cant create an function based on f(x)=ax+b you must simply check if the point A(which you are comparing if its below or above line L)A.xP1x. its above.of if A.x=P1.x then its on it

You could possibly treat it as a graph traversal problem. If you treat adjacent points as 'children' of each other and add each child of your boundary points to a list of open nodes (ie. unvisited nodes), then your enumerator could explore this list.

In case you're unfamiliar with search algorithms, this would probably look something like:
For each point; add the point to a set of closed nodes (ie. visited). If the point is within the polygon, add it to your enumerator and add it's children to your list of open nodes. Continue until your list of open nodes is empty. Also, only add nodes to your open list if they're not already on your closed list.

A couple of notes:

  • I would argue that this approach could be more efficient than simply 'drawing lines' through your shapes. Notably; if any of your polygons are concave, drawing lines may not discover some sections of the shape without checking arbitrarily large numbers of points that aren't within the shape. However, you will end up with a much larger memory overhead than using lines.
  • If you don't need every point straight away, this approach could allow you to generate the next point on-demand.
  • You could improve the performance of the traversal if your problem is amicable to the use of a heuristic. This is most relevant if you're looking for a particular point (or group of points) which satisfy some property.

Have a look at this. It works through some very clear diagrams.

http://www.mathopenref.com/coordpolygonarea2.html

They are just calculating area and you are enumerating points. However you can adapt these ideas. Look carefully at 'A More Complex Case' and how the horizontal lines divide the polygon into a set of trapezoids (if you recognize a triangle as degenerate trapezoid with a zero length side).

If your polygons are (potentially) concave you need to discard some trapezoids based on the direction of travel Y(n)-Y(n+1) rule mentioned in the text.

So, now you've reduced enumerating points in a polygon to enumerating points in a series of trapezoids. That shouldn't be too hard particularly given the trapezoids are nicely oriented being 'trapped' two lines both parallel to the X axis. It is probably easiest to do that in a raster scan so you only have calculate the start X and end X of each raster row.

Your data structure might look like a list of polygons, an index indicating which polygon you're in, which x and y coordinate your at and the end index of the raster row you're in.

I'm assuming you don't mind what order you enumerate in!

Further complexity will be incurred if your 'polygons' included disjoint shapes (i.e. actually more than one polygon or could contain polygonal 'holes'.

The complexity of obtaining the trapezoids is quite low but iterating over a vast number of points may be quite slow. I have no idea what you're doing this for but as pointed out by Stephane Hockenhull above a conventional weight of making it easier is to used a coarsening approximation.

That would mean actually scanning in steps of (say) 5 and (if your algorithm is ameanable) treat the centre point of a 5*5 'cell' as somehow representative.

I hope you weren't looking for an easy answer? I'd love to see some code for this!

Your Answer

 
discard

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.