Skip to content

Conversation

@asinghvi17
Copy link
Member

@asinghvi17 asinghvi17 commented Oct 12, 2025

This still has one Core.Box that I cant seem to get rid of...but is otherwise fully type stable

Benchmarks

Current naive Hao-Sun

Benchmark: 8966 samples with 2 evaluations
 min    13.417 μs
 median 13.604 μs
 mean   14.339 μs
 max    105.875 μs

Prepared with natural tree and using STI

36x speedup!

Benchmark: 8966 samples with 2 evaluations
 min    333.000 ns (2 allocs: 32 bytes)
 median 375.000 ns (2 allocs: 32 bytes)
 mean   403.339 ns (2 allocs: 32 bytes)
 max    18.750 μs (2 allocs: 32 bytes)

Target: TGGeometry

Benchmark: 8966 samples with 2 evaluations
 min    83.000 ns (1 allocs: 16 bytes)
 median 146.000 ns (1 allocs: 16 bytes)
 mean   171.826 ns (1 allocs: 16 bytes)
 max    57.792 μs (1 allocs: 16 bytes)

Benchmarking script

using GADM, NaturalEarth

using GeometryOps: SpatialTreeInterface as STI, 
                   LoopStateMachine,
                   GeometryOps as GO, GeoInterface as GI
import TGGeometry
                
world_fc = naturalearth("admin_0_countries", 110)
world_geoms = GO.tuples(GI.geometry.(GI.getfeature(world_fc)))

world_tree = GO.STRtree(world_geoms)

target_point = (7.5886, 47.5596)

switzerland = GI.geometry(only(GI.getfeature(GADM.get("CHE")))) |> GO.tuples
switzerland_prepared = GO.NaturalIndexing.prepare_naturally(switzerland)

GI.npoint(switzerland)

@assert GO.contains(switzerland, target_point) == true
@assert GO.contains(switzerland, target_point .+ (0, 1)) == false

using Chairmarks

@be GO.contains($(switzerland), $(target_point)), 
    GO.contains($(switzerland_prepared), $(target_point)),
    GO.contains($(GO.TG()), $(GI.convert(TGGeometry, switzerland)), $(GI.convert(TGGeometry, target_point)))

… polygon

This still has one Core.Box that I cant seem to get rid of...
It's cheap and we should be doing it.  Especially because it _really_ helps us going forward.
@asinghvi17
Copy link
Member Author

asinghvi17 commented Oct 19, 2025

So the problem here is the same polygon type issue we ran into earlier. Everything has to be GO.tuplesed before we put it in to any kind of array unless it's pre-sanitized input during polygon clipping.

This sucks but it's the easiest way forward.

p_start = p_end
end

result = SpatialTreeInterface.depth_first_search(per_edge_function,extent -> extent.Y[1] <= y <= extent.Y[2], tree)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The second function could be defined above too this is hard to read

Copy link
Member Author

@asinghvi17 asinghvi17 Oct 19, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah...I might move the internal Hao Sun logic out since it's shared between the tree accelerated and nonaccelerated functions.

if !((v1 < 0 && v2 < 0) || (v1 > 0 && v2 > 0)) # if not cases 11 or 26
u1, u2 = GI.x(p_start) - x, GI.x(p_end) - x
f = Predicates.orient(p_start, p_end, (x, y); exact)
if v2 > 0 && v1 0 # Case 3, 9, 16, 21, 13, or 24
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This wants a little comment giving context to the case numbers

@asinghvi17
Copy link
Member Author

There could be two reasons for the extra overhead:

  • LoopStateMachine overhead (everything returns a LSM.something meaning a bunch of symbol comparisons)
  • Edges are not local to the tree, the C data structure is super optimized for that so the leaf layer of the tree is made up of the edges as an implicit representation, here it's just implied indices.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants