This is mainly to diagnose what is going wrong in your spatial index implementation. Could be useful to add into the package I guess, but keeping it here for now just in case
import GeometryOps as GO
import GeometryOps: SpatialTreeInterface as STI
using GeometryOps.LoopStateMachine: @controlflow
function instrumented_dual_dfs(f::F, predicate::P, node1::T1, node2::T2, log = Any[]) where {F, P, T1, T2}
if STI.isleaf(node1) && STI.isleaf(node2)
# @info "Both leaf"
push!(log, (; n1 = node1, n2 = node2, type = :leaf12, predicate_result = false))
# both nodes are leaves, so we can just iterate over the indices and extents
for (i1, extent1) in STI.child_indices_extents(node1)
for (i2, extent2) in STI.child_indices_extents(node2)
if predicate(extent1, extent2)
@controlflow f(i1, i2)
end
end
end
elseif STI.isleaf(node1) # node2 is not a leaf, node1 is - recurse further into node2
for child in STI.getchild(node2)
if predicate(STI.node_extent(node1), STI.node_extent(child))
push!(log, (; n1 = node1, n2 = child, type = :leaf1, predicate_result = true))
@controlflow instrumented_dual_dfs(f, predicate, node1, child, log)
end
end
elseif STI.isleaf(node2) # node1 is not a leaf, node2 is - recurse further into node1
for child in STI.getchild(node1)
if predicate(STI.node_extent(child), STI.node_extent(node2))
push!(log, (; n1 = child, n2 = node2, type = :leaf2, predicate_result = true))
@controlflow instrumented_dual_dfs(f, predicate, child, node2, log)
end
end
else # neither node is a leaf, recurse into both children
@info "Both not leaf"
for child1 in STI.getchild(node1)
for child2 in STI.getchild(node2)
if predicate(STI.node_extent(child1), STI.node_extent(child2))
push!(log, (; n1 = child1, n2 = child2, type = :inner, predicate_result = true))
@controlflow instrumented_dual_dfs(f, predicate, child1, child2, log)
end
end
end
end
return log
end
This is mainly to diagnose what is going wrong in your spatial index implementation. Could be useful to add into the package I guess, but keeping it here for now just in case