BinnedPointList

mutable struct BinnedPointList{T}

Binned point list structure allowing for fast check for already existing points.

This provides better performance for identifying already inserted points than the naive linear search.

Algorithm

All inserted points are kept in points. In addition, a cuboid binning region (given by binning_region_min/binning_region_max) is subdivided into a number_of_directional_bins-per-dimension grid of bins, each holding the indices of the points falling into it. Looking up or inserting a point then only has to search the (small) bin the point falls into, instead of the whole point list.

Points that fall outside of the current binning region (or are inserted before any region has been established) are kept in the separate unbinned list, which is searched linearly. Once too many points accumulate there (as governed by num_allowed_unbinned_points and max_unbinned_ratio), the binning region is enlarged to cover them and all points are re-distributed into freshly allocated bins, see _rebin_all_points!. This amortizes the cost of handling points outside the initial region over many insertions.

OTOH the implementation is still quite naive - it dynamically maintains a cuboid binning region with a fixed number of bins.

Probably tree based adaptive methods (a la octree) will be more efficient, however they will be harder to implement.

In an ideal world, we would maintain a dynamic Delaunay triangulation, which at once could be the starting point of mesh generation which will follow here anyway.

• dim::Int32: Space dimension

• tol::Any: Point distance tolerance. Points closer than tol (in Euclidean distance) will be identified, i.e. are collapsed to the first inserted. This is an absolute, not a relative, tolerance.

• binning_region_min::Vector: The union of all bins is the binning region - a cuboid given by two of its corners. It is calculated dynamically depending on the inserted points. Until the first call to _rebin_all_points! actually performs a rebinning, this is the (degenerate, empty) cuboid binning_region_min = fill(floatmax(T), dim), binning_region_max = fill(-floatmax(T), dim).

• binning_region_max::Vector

• binning_region_increase_factor::Any: Increase factor of binning region (with respect to the cuboid defined by the coordinates of the binned points). Keeps points which lie exactly on the boundary of the region spanned by the current point set (and future points close to it) safely inside the binning region instead of the unbinned list.

• points::ElasticArrays.ElasticArray{T, 2, M, V} where {T, M, V<:DenseVector{T}}: The actual point list, as a dim × npoints matrix.

• bins::Array{Vector{Int32}}: The bins are vectors of indices of points in the point list. We store them in a dim-dimensional array of size number_of_directional_bins along each dimension, i.e. with number_of_directional_bins^dim entries in total. Before the first rebinning, this array has size zero in each dimension (no bins have been allocated yet).

• number_of_directional_bins::Int32: Number of bins in each space dimension

• unbinned::Vector{Int32}: Some points will fall outside of the binning region, e.g. because they were inserted before the region was ever calculated, or because they lie outside of the region as it was last enlarged. We collect them in a vector of unbinned point indices, which is searched linearly.

• num_allowed_unbinned_points::Int32: Number of unbinned points tolerated without triggering a rebinning via _rebin_all_points!.

• max_unbinned_ratio::Any: Maximum ratio of unbinned points in the point list tolerated without triggering a rebinning via _rebin_all_points!. Rebinning is triggered once the number of unbinned points exceeds max(num_allowed_unbinned_points, max_unbinned_ratio * npoints).

• current_bin::Vector{Int32}: Storage of current point bin

findpoint(binnedpointlist, p)

Find point p in binned point list. Return its index in the point list if found, otherwise return 0. p may be a vector or a tuple.

This first triggers a rebinning if too many points have accumulated in the unbinned list (see _rebin_all_points!), then looks up the bin p falls into (or the unbinned list, if p falls outside of the current binning region) and searches it linearly.

 Base.insert!(binnedpointlist,p)

If another point with distance less than tol from p is already in the point list, return its index. Otherwise, insert p into the point list and return the index of the newly inserted point. p may be a vector or a tuple.

 Base.insert!(binnedpointlist,x)

Insert 1D point via coordinate.

 Base.insert!(binnedpointlist,x,y,z)

Insert 3D point via coordinates.

_findpoint(binnedpointlist, index, p)

Find point p by linear search among the points listed in index (typically one bin, or the unbinned list). Return its index in binnedpointlist.points, or zero if not found.

_bin_of_point!(binnedpointlist, p)

Calculate the bin of the point. Result is stored in bpl.current_bin.

For each space dimension, bpl.current_bin[idim] is set to a value in 1:bpl.number_of_directional_bins if p[idim] lies within the binning region for that dimension, and to 0 if it lies outside (this includes the case where the binning region has not been set up yet, or is still degenerate in that dimension). Consumers must check reduce(*, bpl.current_bin) > 0 before indexing into bpl.bins with bpl.current_bin, since a 0 entry is not a valid array index.

_rebin_all_points!(bpl)

Re-calculate binning if there are too many unbinned points, i.e. if length(bpl.unbinned) > max(bpl.num_allowed_unbinned_points, bpl.max_unbinned_ratio * size(bpl.points, 2)). This amounts to two steps:

• Enlarge the binning region in order to include all (previously unbinned) points, plus some tolerance-dependent slack given by binning_region_increase_factor so that points exactly on the boundary are unambiguously inside.
• Re-allocate bpl.bins and re-calculate the bin of every point (not just the previously unbinned ones), since the enlarged region shifts the boundaries of all bins.

This is a full O(npoints) operation; it is only called when the above threshold is exceeded, so its cost is amortized over the insertions/lookups performed since the last rebinning.

naiveinsert!(binnedpointlist, p)

Insert via linear search, without any binning. Just for being able to check if all of the above was worth the effort...