Linked Lists Grid Structure for fixed-radius neighbor search

After reading Optimization of Algorithm to Fill Empty Space, I decided to try and implement the proposed algorithm (see this response) for fixed-radius neighbor search.

The goal is, given a set of points in 2D or 3D space, to find all pairs of points whose euclidean distance are inferior to a given radius.

The proposed idea is to use a uniform grid structure:

  • The space is partitioned into cells (squares in 2D, cubes in 3D) of equal size, whose width correspond to the search radius
  • Each point is assigned to its cell
  • Now, to find neighbors of a given point, one only needs to search within the cells that are adjacent to the cell of the point.

I followed the implementation proposed here, based on linked lists. The grid is represented as

  • a head of type Vec<Option<usize>>, whose length is the number of cells, containing the index of the first element of each cell.
  • a successor of type Vec<Option<usize>>, whose length is the number of elements in the grid, storing the successor of each element in the cell.

Now here's my rust implementation for a 3D space. I commented it as much as possible, so hopefully the code should be clear if you understood the implementation details above. (if not, let me know !)

I'd love to get feedback since I am pretty new to rust. In particular:

  • Was creating new structures to iterate over cells and neighbors the right approach ?
  • Is this linked list structure optimal in terms of performance ? If not, what would be a better one ?
  • Is there a way I can make this liked list implementation more performant ?
  • How should I go about processing / generating the iterator of neighbors in parallel ?
use crate::itertools::{Itertools, Product};

use smallvec::SmallVec;
use std::iter::Iterator;

use nalgebra::{distance_squared, Point3};

const DIM: usize = 3;

type CellIdx = usize;

type EntityIdx = usize;

pub struct GridPartition {
    // squared search radius
    squared_radius: f32,
    // number of subdivisions for each dimension
    pub nb_subdivisons: [usize; DIM],
    // length of the cells for each dimension (slightly larger than the search radius)
    pub cell_width: [f32; DIM],
    // total number of cells
    pub nb_cells: usize,
    // adjascent_cells[cell_idx] = indexes of all cells adjascent to 'cell_idx'
    // adjascent_cells.len() = #cells
    pub adjascent_cells: Vec<SmallVec<[CellIdx; 14]>>,
    // head[cell_idx] = position of the first element of the cell in self.successor
    // head.len() = #cells.
    head: Vec<Option<EntityIdx>>,
    // successor[entity_idx] = idx of the sucessor of entity 'entity_idx' in the cell
    // successor.len() = #entities.
    successor: Vec<Option<EntityIdx>>,

impl GridPartition {
    // Create a new grid, provided a search radius
    // Note: the points are assumed o belong to 
    // [0, dimensions[0]] x [0, dimensions[1]] x [0, dimensions[2]] 
    pub fn new(dimensions: [f32; DIM], radius: f32) -> Self {
        let squared_radius = radius.powi(2);
        let nb_subdivisons = [
            (dimensions[0] / radius) as usize,
            (dimensions[1] / radius) as usize,
            (dimensions[2] / radius) as usize,
        let cell_width = [
            dimensions[0] / nb_subdivisons[0] as f32,
            dimensions[1] / nb_subdivisons[1] as f32,
            dimensions[2] / nb_subdivisons[2] as f32,
        let nb_cells = (nb_subdivisons[0] * nb_subdivisons[1] * nb_subdivisons[2]) as usize;
        let adjascent_cells = vec![SmallVec::<[CellIdx; 14]>::new(); nb_cells];
        let head = vec![None; nb_cells];
        let successor = Vec::new();
        let mut grid = Self {

    // Precompute the 'adjascent_cells' vector
    fn build_adjascent_cells(&mut self) {
        // for every cell
        for idx_flat in 0..self.nb_cells {
            let idx_spatial = self.unflatten_index(idx_flat);
            // for all adjascent cells
            for &dx in &[-1, 0, 1] {
                for &dy in &[-1, 0, 1] {
                    for &dz in &[-1, 0, 1] {
                        // if the adjascent cell is not out of bounds
                        if (idx_spatial.x as i64 + dx >= 0)
                            && (idx_spatial.x as i64 + dx < self.nb_subdivisons[0] as i64)
                            && (idx_spatial.y as i64 + dy >= 0)
                            && (idx_spatial.y as i64 + dy < self.nb_subdivisons[1] as i64)
                            && (idx_spatial.z as i64 + dz >= 0)
                            && (idx_spatial.z as i64 + dz < self.nb_subdivisons[2] as i64)
                            let other_idx = Point3::new(
                                (idx_spatial.x as i64 + dx) as usize,
                                (idx_spatial.y as i64 + dy) as usize,
                                (idx_spatial.z as i64 + dz) as usize,
                            let other_idx_flat = self.flatten_index(other_idx);
                            // makes sure to avoid duplicate neighbors: two cells (i, j)
                            // are adjascent only if i <=j
                            if other_idx_flat <= idx_flat {
                                // add the adjascent cell to the vector

    // given the sptial coordinates (i, j, k)  of a cell, compute its integer index.
    fn flatten_index(&self, idx_spatial: Point3<usize>) -> usize {
        idx_spatial.x * self.nb_subdivisons[1] * self.nb_subdivisons[2]
            + idx_spatial.y * self.nb_subdivisons[2]
            + idx_spatial.z

    // given the integer index of a cell, compute its (i, j, k) coordinates
    fn unflatten_index(&self, idx_flat: CellIdx) -> Point3<usize> {
            idx_flat / (self.nb_subdivisons[1] * self.nb_subdivisons[2]),
            (idx_flat / self.nb_subdivisons[2]) % self.nb_subdivisons[1],
            idx_flat % self.nb_subdivisons[2],

    // Return an iterator over the entities contained in a cell
    pub fn iter_cell(&self, idx_flat: CellIdx) -> CellIterator {
        CellIterator {
            current_entity: self.head[idx_flat],
            successor: &self.successor,

    // Return an iterator over the entities contained in cells adjascent to cell 'idx_flat'
    pub fn iter_adjascent_cells(&self, idx_flat: CellIdx) -> MultiCellIterator {
        let adjascent_cells = &self.adjascent_cells[idx_flat];
        let current_cell = Some(self.iter_cell(adjascent_cells[0]));
        let remaining_cells: SmallVec<[CellIterator; 14]> = adjascent_cells[1..]
            .map(|&i| self.iter_cell(i))
        MultiCellIterator {

    // Return an iterator over all couples the couples (i, j), where j is in cell 'idx_flat',
    // and i is in an cell adjascent to 'idx_cell'
    fn iter_cell_couples(&self, idx_flat: CellIdx) -> Product<MultiCellIterator, CellIterator> {

    // Return an iterator over all couples (i, j), where i and j are within the (squared) distance self.squared_radius
    pub fn iter_neighbors<'a>(&'a self, points: &'a [Point3<f32>]) -> NeighborIterator {
        NeighborIterator {
            grid: self,
            current_cell: 0,
            cartesian_iterator: self.iter_cell_couples(0),

    // Reset the head and the successor vector
    fn empty_grid(&mut self, nb_points: usize) {
        // Reset the 'head' vector
        for x in &mut self.head {
            *x = None;
        // Resize the 'successor' vector if needed
        let additional_capacity = nb_points - self.successor.capacity();
        if additional_capacity > 0 {
        // Fill the successor vector with None values.
        // the values already present are not replaced,
        // since they will be overwritten in self.fill
        let additional_points = nb_points - self.successor.len();
        if additional_points > 0 {
            for _ in 0..additional_points {

    // fill the grid with a new set of points.
    fn fill_grid(&mut self, points: &[Point3<f32>]) {
        // for each point
        for (i, point) in points.iter().enumerate() {
            let idx_spatial = Point3::new(
                (point.x / self.cell_width[0]) as usize,
                (point.y / self.cell_width[1]) as usize,
                (point.z / self.cell_width[2]) as usize,
            let idx_flat = self.flatten_index(idx_spatial);
            // define the successor of point as the previous first element of the cell
            self.successor[i] = self.head[idx_flat];
            // define the first element of the cell as the current point
            self.head[idx_flat] = Some(i);

    // Update the grid with a new set of points,
    // and return an iterator over all the neighbor indexes
    pub fn query_neighbors<'a>(&'a mut self, points: &'a [Point3<f32>]) -> NeighborIterator {

// Iterator over all the entities in a cell
pub struct CellIterator<'a> {
    current_entity: Option<EntityIdx>,
    successor: &'a Vec<Option<EntityIdx>>,

impl<'a> Iterator for CellIterator<'a> {
    type Item = EntityIdx;

    fn next(&mut self) -> Option<EntityIdx> {
        // if their is a current entity
        if let Some(i) = self.current_entity {
            // get the next entity
            self.current_entity = self.successor[i];
            // return the current entity
        // else: we finished iterating the cell
        } else {

// Iterator over the entities of several cells
pub struct MultiCellIterator<'a> {
    // Would it be better to use an iterator for 'remaining_cells' instead of an array ? How to do it ?
    current_cell: Option<CellIterator<'a>>,
    remaining_cells: SmallVec<[CellIterator<'a>; 14]>,

impl<'a> Iterator for MultiCellIterator<'a> {
    type Item = EntityIdx;

    fn next(&mut self) -> Option<EntityIdx> {
        // if their is a current cell to iterate over
        if let Some(cell) = &mut self.current_cell {
            // if the current cell is not empty
            if let Some(i) = &mut {
                // pop the element of the cell
            // else: we finished iterating the current cell
            } else {
                // go to the next cell
                self.current_cell = self.remaining_cells.pop();
                // keep iterating
        // else: we have finished visiting all cells
        } else {

// Iterator over all neighbors of the Grid
pub struct NeighborIterator<'a> {
    grid: &'a GridPartition,
    current_cell: CellIdx,
    cartesian_iterator: itertools::Product<MultiCellIterator<'a>, CellIterator<'a>>,
    points: &'a [Point3<f32>],

impl Iterator for NeighborIterator<'_> {
    type Item = (CellIdx, EntityIdx, f32);

    fn next(&mut self) -> Option<(EntityIdx, EntityIdx, f32)> {
        // if the cartesin iterator is not empty
        if let Some((idx, idx_other)) = {
            if idx < idx_other {
                // if the points are close enought, yield the neighbors and their distance
                let squared_distance =
                    distance_squared(&self.points[idx], &self.points[idx_other]) as f32;
                if squared_distance < self.grid.squared_radius {
                    Some((idx, idx_other, squared_distance.powf(0.5)))
                } else {
            // otherwise, keep iterating
            } else {
        // else: the cartesian iterator is empty
        } else {
            // go to the next cell
            self.current_cell += 1;
            // if the current cell is not the last one
            if self.current_cell < self.grid.nb_cells {
                // compute the cartesian iterator of the next cell
                self.cartesian_iterator = self.grid.iter_cell_couples(self.current_cell);
                // keep iterating
            // else:
            } else {
                // the cartesian iterator is empty this is the last cell
                // we have finised iterating

Thank you if you took the time to read through the code :slight_smile:

1 Like

To keep the thread organized, I have chosen to comment here.

Thank you for sharing your implementation. I am looking at the code right now and have looked for a bit, but I am struggling with;

  1. How do I initialize a grid?
  2. How do I fill the grid with neighbours?
  3. How do I extract particles in a neighbour?

The only thing I have got to work right now is:

const DIM: usize = 3;
let grid = celllist::GridPartition::new([1.0; DIM],0.08);

But I don't know how to make it insert the points etc.

I really like your implementation and would like to help improve it, but I need a bit of a helping hand with you making a simple example of some random points scattered in space if possible.

Once again thank you very much, I would really like to use it, but need a bit of help

Kind regards

Just out of curiosity, what sorts of applications is this fixed-size search typically used for? I do a lot of work with CAD software and whenever I need to do spatial-based lookups it'll almost always be done based on user-provided input.

I remember watching a talk by Bjarne Stroustrup on the performance of linked lists versus good ol' std::vector<T>. There's also a nice explanation on StackOverflow.

They showed how linked lists are a terrible data structure to use if you want to do any sort of iteration, and even for operations that linked lists are algorithmically faster for (e.g. inserting new elements in the middle), the overhead of doing all that pointer chasing and cache missing completely dominates any algorithmic improvements you get.

Also, processors are ludicrously fast at doing linear scans (e.g. Vec<T>)... when it sees your code accessing memory sequentially the prefetcher will automatically load that memory into cache in the background so that every lookup comes from cache.

It also decreases the memory usage because each item doesn't have the overhead of one (or two for doubly-linked lists) extra pointers you get with linked lists. This can sometimes affect performance because you are able to fit more items into cache.

TL;DR: I'd just use a Vec<T> instead of a linked list.

In general, if your type is Sync and has iterator types for the operations you're wanting to do, you're good to go. Just try it, and if the compiler doesn't yell at you then it will all Just Work. Such is the power of the borrow checker :slightly_smiling_face:

You probably won't be able to parallelize things like fill_grid() which mutate the grid, though. To do that you'd need a way to make sure each thread only populates its own section of the grid, and that no sections overlap. I feel like this would be tricky to do in a way that the borrow checker is satisfied and doesn't require unsafe code to manage borrowing at runtime.

Something else to keep in mind is that parallelism works best when done coarsely. For example, parallelism would be great have different "tasks" where you need to do a bunch of searches and can run each task on its own thread. But when the operation each thread is doing is small (e.g. incrementing a number) then the overhead of multithreading will start to dominate.

Of course, the best advice when doing any sort of performance tuning is to write realistic benchmarks and use them to compare the effect of different algorithms/data structures/memory layout. Consider the above a rule of thumb, there's no guarantee it'll be valid for your use case.

Here is an example of how to use it.

The method query builds the grid (through fill_grid) given a slice of points, and returns an iterator over the indexes of neighbors as well as their distance (through iter_neighbors).

use rand::Rng;

fn random_points(n: usize) -> Vec<Point3<f32>> {
    let mut rng = rand::thread_rng();
    let mut x = Vec::with_capacity(n);
    for _ in 0..n {
            rng.gen_range(0., 1.),
            rng.gen_range(0., 1.),
            rng.gen_range(0., 1.),

fn main() {
    let dimensions = [1., 1., 1.];
    let radius = 0.1;
    let nb_points = 100;
    let points = random_points(nb_points);
    let mut grid = GridPartition::new(dimensions, radius);
    for (i, j, d) in grid.query_neighbors(&points) {
            "{:?} is close to {:?} => distance = {}",
            points[i], points[j], d

This fixed size search is typically used for simulations involving many particles / entities, coupled by a short range interactions, such as particle-based fluid simulations.

I saw warnings about linked lists in rust's standard library so I was more or less aware of the issue you mentioned. I am currently implementing a vector-based version of the algorithm, however it requires sorting the array of entities according to the cell they belong to (making the overall algorithm complexity nlog(n) rather than n). I'll post a comparison here once it's done.

The vector-based implementation is explained here (Dynamic Hashed Grid section) in case anyone is interested !

Regarding parallelization, creating the grid is quite fast and isn't the performance bottleneck, so it shouldn't be too hard to do. I'll look into how to do that using Sync, thanks :slight_smile:


I will try to keep playing around with it until I grasp it even better. Didn't know that nalgebra had a Point class already defined so will try to change my old code to accomodate that too.

Also very interesting with the hashing etc. using vector based implementation. Hoping you are able to achieve it, I know it would benefit me personally a lot, since I have always struggled with these neighbour search things.

Kind regards

This topic was automatically closed 90 days after the last reply. New replies are no longer allowed.