How to pre-process 2D polygon to speed up checking if a point is inside a polygon

[DuoZhangRobotics]

Issue Details

Hi,

Thanks in advance for all your help. My code is currently checking whether a point is inside a polygon or not intensively on the same set of polygons. I am wondering if there are any methods that allow me to preprocess the polygon beforehand so I can do the checking a lot faster later. Thanks again.

Duo

Environment

• Operating system (Windows/Mac/Linux, 32/64 bits): Ubuntu 22.04
• Compiler: gcc11.4
• Release or debug mode: Release
• Specific flags used (if any):
• CGAL version: 5.6.1
• Boost version: 1.74
• Other libraries versions if used (Eigen, TBB, etc.):

[afabri]

You can create a 2D constrained triangulation, and mark the triangles that are inside the polygons. See this example. You then call the cdt.locate(p) function. If you know all points do a spatial_sort() with them, and when you do locate() give the face handle returned by the previous locate() as hint.

[DuoZhangRobotics]

Thanks! I'll try it out!

[afabri]

I did not understand what you mean with "or not intensively on the same set of polygons." Maybe a drawing would help.

[DuoZhangRobotics]

Hi Andreas,

I have a bunch of polygons that are already fixed. And I have a huge set of points. I want to know whether a point is inside of a polygon or not, for every pair of point and polygon. Since all the polygons will be queried multiple times, I was wondering if there's any method to pre-process the polygons so I can speed up all the checking.

[efifogel]

Insert all the polygons into a 2D arrangement data structure, and then apply point-location queries. You can even apply batch point location operations, where you use a bunch of points and you get the location of each in an efficient way.

---

/_____/) o /_________ __ // (____ ( ( ( (/ (/-(-'_(/ _/

On Sat, 19 Oct 2024 at 00:15, Duo Zhang @.***> wrote:

Hi Andreas,

I have a bunch of polygons that are already fixed. And I have a huge set of points. I want to know whether a point is inside of a polygon or not, for every pair of point and polygon. Since all the polygons will be queried multiple times, I was wondering if there's any method to pre-process the polygons so I can speed up all the checking.

— Reply to this email directly, view it on GitHub https://github.com/CGAL/cgal/issues/8549#issuecomment-2423242975, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABVBNODFEJYM5W3XATWTLXTZ4F26ZAVCNFSM6AAAAABP7QYQBKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDIMRTGI2DEOJXGU . You are receiving this because you are subscribed to this thread.Message ID: @.***>

[raphaelsulzer]

If I understand correctly, the OP does not only want to know if a point is inside any polygon, but also inside which one. Say, for example, you want to assign LiDAR points to a set of 2D building footprint polygons.

If I am not mistaken, it is not possible to retrieve the input polygons from face handles returned by the locate function of the 2D constraint triangulation / 2D arrangement?

Something like an AABB tree of a set of polygons could help here. But that does not exist either, right?

[efifogel]

Of course you can. It takes a bit more care and processing. Extend the face record of the underlying halfedge data structure with an identifier of the polygon.

---

/_____/) o /_________ __ // (____ ( ( ( (/ (/-(-'_(/ _/

On Sat, 19 Oct 2024 at 13:09, Raphael Sulzer @.***> wrote:

If I understand correctly, the OP does not only want to know if a point is inside any polygon, but also inside which one. Say, for example, you want to assign LiDAR points to a set of 2D building footprint polygons.

If I am not mistaken, it is not possible to retrieve the input polygons from face handles returned by the locate function of the 2D constraint triangulation / 2D arrangement.

Something like an AABB tree of a set of polygons could help here. But that does not exist either, right?

— Reply to this email directly, view it on GitHub https://github.com/CGAL/cgal/issues/8549#issuecomment-2423739697, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABVBNOEH6EUAXAZBNQIIZ7TZ4IVXBAVCNFSM6AAAAABP7QYQBKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDIMRTG4ZTSNRZG4 . You are receiving this because you commented.Message ID: @.***>

[afabri]

I guess it needs more input from @DuoZhangRobotics Are the polygons small but spread over a big area? Do all polygons intersect many other polygons? In case there is a representative data set it might be interesting if you shared it with us.

[afabri]

[DuoZhangRobotics]

Sorry for the late reply. These are two simple sets of polygons in my case. They are basically some visible areas from the CGAL's 2D Visiblity library. If I use the locate functions from CGAL, it seems a little bit too slow. If I convert the arrangements of the polygons into the polygon class in CGAL, and then call cgal::bounded_side_2() function, it'll be so much faster. But in this way, I can just use plain C++ implementations with float type, then I guess it'll be faster than CGAL in this case. And I don't think I can convert all the polygons to an arrangement and do batch checking here? Because the polygons have too many intersections among them. In this case, AABB tree or BVH wouldn't help too much, but there will be some improvements. I'm also thinking about using Cuda or avx to parallelize everything. But according to my benchmarking, the AVX code doesn't help much when -O3 is turned on compared with my plain c++ code, but they are both approximately 2times faster than bounded_side_2() function. I guess that's because I'm using the exact kernel in CGAL and only using float for AVX and C++ only code.

[afabri]

I have a hard time to single out the different polygons. Can you give them different colors and maybe not more than 10 just for identifying them.

[efifogel]

On Sun, 20 Oct 2024 at 18:07, Duo Zhang @.***> wrote:

image.png (view on web) https://github.com/user-attachments/assets/80904cab-9485-4279-983c-9e747969c62e image.png (view on web) https://github.com/user-attachments/assets/9bbf8f60-57b0-4fe4-b8ae-cd0298f18268 Sorry for the late reply. These are two simple sets of polygons in my case. They are basically some visible areas from the CGAL's 2D Visiblity library. If I use the locate functions from CGAL, it seems a little bit too slow. If I convert the arrangements of the polygons into the polygon class in CGAL, and then call cgal::bounded_side_2() function, it'll be so much faster.

No it will not, or better stated, point-location in 2D arrangements is extremely fast. As a matter of fact, 2D arrangements are used as the underlying data structure of most Boolean operations on polygons.

But in this way, I can just use plain C++ implementations with float type, then I guess it'll be faster than CGAL in this case. And I don't think I can convert all the polygons to an arrangement and do batch checking here? Because the polygons have too many intersections among them.

If you can afford spending time on preprocessing, which consists of constructing the arrangement, then you will benefit from the efficient point-location. Construction is done with a sweep line alg,. which is output sensitive and depends on the number of intersections (O(n log k), where n is the number of edges of all polygons and k is the number of intersections in your case.)

In this case, AABB tree or BVH wouldn't help too much, but there will be

some improvements. I'm also thinking about using Cuda or avx to parallelize everything. But according to my benchmarking, the AVX code doesn't help much when -O3 is turned on compared with my plain c++ code, but they are both approximately 2times faster than bounded_side_2() function. I guess that's because I'm using the exact kernel in CGAL and only using float for AVX and C++ only code.

— Reply to this email directly, view it on GitHub https://github.com/CGAL/cgal/issues/8549#issuecomment-2425034755, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABVBNOHBXM7QBSK6RJEIAO3Z4PBMVAVCNFSM6AAAAABP7QYQBKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDIMRVGAZTINZVGU . You are receiving this because you commented.Message ID: @.***>

[efifogel]

Consider splitting the input set of polygons into subsets of polygons, such that for each subset the number of intersections is small. (This may not be easy or even possible...)

---

/_____/) o /_________ __ // (____ ( ( ( (/ (/-(-'_(/ _/

On Sun, 20 Oct 2024 at 18:28, Efi Fogel @.***> wrote:

On Sun, 20 Oct 2024 at 18:07, Duo Zhang @.***> wrote:

image.png (view on web) https://github.com/user-attachments/assets/80904cab-9485-4279-983c-9e747969c62e image.png (view on web) https://github.com/user-attachments/assets/9bbf8f60-57b0-4fe4-b8ae-cd0298f18268 Sorry for the late reply. These are two simple sets of polygons in my case. They are basically some visible areas from the CGAL's 2D Visiblity library. If I use the locate functions from CGAL, it seems a little bit too slow. If I convert the arrangements of the polygons into the polygon class in CGAL, and then call cgal::bounded_side_2() function, it'll be so much faster.

No it will not, or better stated, point-location in 2D arrangements is extremely fast. As a matter of fact, 2D arrangements are used as the underlying data structure of most Boolean operations on polygons.

But in this way, I can just use plain C++ implementations with float type, then I guess it'll be faster than CGAL in this case. And I don't think I can convert all the polygons to an arrangement and do batch checking here? Because the polygons have too many intersections among them.

If you can afford spending time on preprocessing, which consists of constructing the arrangement, then you will benefit from the efficient point-location. Construction is done with a sweep line alg,. which is output sensitive and depends on the number of intersections (O(n log k), where n is the number of edges of all polygons and k is the number of intersections in your case.)

In this case, AABB tree or BVH wouldn't help too much, but there will be

some improvements. I'm also thinking about using Cuda or avx to parallelize everything. But according to my benchmarking, the AVX code doesn't help much when -O3 is turned on compared with my plain c++ code, but they are both approximately 2times faster than bounded_side_2() function. I guess that's because I'm using the exact kernel in CGAL and only using float for AVX and C++ only code.

— Reply to this email directly, view it on GitHub https://github.com/CGAL/cgal/issues/8549#issuecomment-2425034755, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABVBNOHBXM7QBSK6RJEIAO3Z4PBMVAVCNFSM6AAAAABP7QYQBKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDIMRVGAZTINZVGU . You are receiving this because you commented.Message ID: @.***>

[DuoZhangRobotics]

Hi @afabri, sorry for the late reply again. I was kinda swamped by some errands for last two days. Here's another screenshot with only 6 polygons with different colors. Hope this will help.

Hi @efifogel, thanks for your input too! If I build all the polygons into one arrangement, given a new point, I should still call bounded_side_2 right? If so, then can I still tell which polygon the point lands in?

I really appreciate all the help!

[raphaelsulzer]

I think what @efifogel is suggesting is

1. Insert your polygons into the arrangement (see e.g. point 3. here). The insert functions return Halfedge handles that you associate with your polygons. Either through a simple map, or by extending the Halfedge type of the Arrangement.
2. Then, do point location queries with your query points. Most of the points will be located in a face. From the Face handle you can iterate over the face's halfedges. Finally, from the halfedges you can look up/know in which input polygon you are in.

This workflow assumes your input polygons are not overlapping. If they are, it is still possible, but needs some intermediate steps.

[efifogel]

This is exactly what I meant. There are several ways of inserting polygons into an arrangement and maintaining an association between the cells and additional data, such as labels. Attached is a simple method, which may not be the most efficient. The construction of the arrangement is the dominant part, so if you care about performance, you may want to implement a different method.

1. In the implemented method the polygons are assumed to be open. If you need them closed, you will also need to extend the halfedge and vertex records and implement all the 9 nops operations in the overlay traits.

2. Another method is inserting the segments that comprise the polygon boundaries directly into the final arrangement and use an observer to keep track of the labels of the halfedges. Then you need to traverse all faces in a BFS or a DFS manner. When you visit a new face you add the label(s) obtained from the haldge you just crossed, to the current list of labels, and use them to update the label(s) in the newly visited face. When you exit a visited face you remove the label(s) obtained from the haldge you just crossed.

   ---

/_____/) o /_________ __ // (____ ( ( ( (/ (/-(-'_(/ _/

On Fri, 25 Oct 2024 at 01:16, Raphael Sulzer @.***> wrote:

I think what @efifogel https://github.com/efifogel is suggesting is

1. Insert your polygons into the arrangement (see e.g. point 3. here https://doc.cgal.org/latest/Arrangement_on_surface_2/index.html#aos_sec-tips. The insert functions return Halfedge handles that you associate with your polygons. Either through a simple map, or by extending the Halfedge type of the Arrangement.
2. Then, do point location queries https://doc.cgal.org/latest/Arrangement_on_surface_2/index.html#title17 with your query points. Most of the points will be located in a face. From the Face handle you can iterate over the face's halfedges. Finally, from the halfedges you can know in which input polygon you are in. This workflow assumes your input polygons are not overlapping. If they are, it is still possible, but needs some intermediate steps.

— Reply to this email directly, view it on GitHub https://github.com/CGAL/cgal/issues/8549#issuecomment-2436435577, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABVBNOAVNO5XQYKO2M6JCLLZ5FWS3AVCNFSM6AAAAABP7QYQBKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDIMZWGQZTKNJXG4 . You are receiving this because you were mentioned.Message ID: @.***>

cmake_minimum_required(VERSION 3.29.0) project(ppl) find_package(CGAL QUIET COMPONENTS Core OPTIONAL_COMPONENTS Qt6) find_package(Boost REQUIRED COMPONENTS filesystem program_options system regex)

if (NOT CGAL_FOUND) message(FATAL_ERROR "This program requires the CGAL library, and will not be compiled.") endif()

if (CMAKE_COMPILER_IS_GNUCXX) add_definitions(-std=c++17) endif()

create_single_source_cgal_program("ppl.cpp") target_link_libraries(ppl PUBLIC ${Boost_LIBRARIES})

if(CGAL_Qt6_FOUND) target_link_libraries(ppl PUBLIC CGAL::CGAL_Basic_viewer) endif()

[DuoZhangRobotics]

Thanks so much! I think I can give it a try. Thanks for all your help!