Guide

What is worth knowing before the reference pages. Everything here applies to both APIs unless it says otherwise.

Indices are stable, and errors are exceptions

Three things to know before using it. A shell is Rmin <= d <= Rmax, closed at both ends, so a ball holds every object within Rmax and an object sitting exactly on either radius is returned. Indices are stable: removeObject leaves an index permanently invalid, addObject issues a new one, and indices are never reused or renumbered, so a stored index never comes to mean a different object. And errors are exceptions — meshsearch::Error, or meshsearch::IndexError for an index that is out of range or removed — never return codes.

The overloads taking an object exclude it by identity, not by distance: closeObjects(i, Rmax) omits object i itself, but another object at the same position is a different object, is at d = 0, and is returned. The same holds for nearestObject and nearestObjects. The point-taking overloads have no identity to exclude, so an object lying on the query point is simply within the radius and is returned.

The constructor checks its input before allocating anything. It raises meshsearch::Error for a cellsize that is not finite and positive, for limits that are not finite or have lo > hi, for a coordinate that is not finite, and for a grid too large to count: more than 2147483647 cells along an axis, or more than 4294967295 objects. Running out of memory is not translated: it is std::bad_alloc in C++ and MemoryError in Python, and an addObject that fails that way leaves the grid as it was.

Const members may be called concurrently on one grid; non-const members may not.

The cell side

The cell side is the one tuning parameter. Too large and every query tests the many objects sharing a cell; too small and a query walks many cells to cover the same sphere. One to four mean separations of the data is the usual range, and two to three is a reasonable default for mixed use.

What each choice costs is measured, at two data sizes and for all three query kinds, on the performance page. For a shell query the number that decides the work is Rmax / cellsize rather than either alone: keep it at or below 1.

The constructor takes the box as well, and if it is not given, fits it to the data and grows each side to a whole number of cells — so get_lims (lims in Python) can report a box slightly wider than the points, and a point just outside the data range may still be accepted.

An axis along which the box has no extent — every point at the same coordinate, or limits with lo == hi — gets one cell, and the box is grown to one cell side around it. Flat and needle-shaped point sets are therefore accepted. A box that cannot be counted is refused with an error instead: more than 2147483647 cells along one axis, which a cell side tiny against the box would ask for.

Batched queries, and the form they return

The three point-taking queries also accept arrays of query points, answering all of them in one call:

nearest = grid.nearest_object(qx, qy, qz)              # (n_query,) uint32
five    = grid.nearest_objects(5, qx, qy, qz)          # (n_query, 5) uint32
idx, off = grid.close_objects(qx, qy, qz, 0.02)        # compressed-row form

The same names are overloaded: three scalars ask about one point, three arrays ask about many. Query points must be numpy arrays rather than lists, because an array argument is what selects the batched form. A one-element array takes the batched path and returns a one-element result.

This exists only in Python, and that is deliberate. Crossing the Python boundary costs a few microseconds per call, which is more than some of these queries take to run: it is why calling scipy’s cKDTree once per point is about five times slower than handing it the whole array, on identical work. A C++ caller has no such boundary and writes a loop, so the C++ library has no batched API and is not missing one. The index-taking overloads have no batched form either, in Python or C++: their input is a single index, so there is no per-call cost to amortise over anything.

A batched call raises where the single-point call would, naming the query that failed — query 17: X coordinate outside the box — and returns nothing partial.

The compressed-row result

nearest_objects is rectangular: every query returns exactly n, so the result is one (n_query, n) array. A shell is not — the count varies per query — so close_objects returns two arrays instead:

• indices holds every result, concatenated in query order.

• offsets has n_query + 1 entries. The results of query i are indices[offsets[i]:offsets[i+1]]. offsets[0] is 0 and the last entry is len(indices). A query with no results gives an empty slice, never a missing entry. Both arrays are uint32, and offsets counts results, not queries, so it is wide enough for the whole batch.

idx, off = grid.close_objects(qx, qy, qz, 0.02)

idx[off[i]:off[i+1]]          # the neighbours of query i, no copy
counts = np.diff(off)         # how many each query found
counts.sum() == len(idx)      # always

The idiom that makes this form worth learning is pairing each result with the query it came from:

query_of = np.repeat(np.arange(len(qx)), np.diff(off))

indices and query_of are now two flat arrays of the same length, so any per-query quantity is one vectorised expression over every result at once — np.bincount(query_of, weights=...) to reduce per query, boolean masks to select across queries. examples/example_batch.py works through it.

Why not a list of arrays? Because that would allocate one Python object per query, which is the cost the batched call exists to remove: a batch of 50,000 queries would hand back 50,000 numpy arrays. There is deliberately no helper to convert the compressed form into such a list, it would be slow and would undo the point. Slice it, or use the flat form directly.