I am dealing with a set of $N$ dimensionless (point) particles in a box. The box has a certain volume $V$. I need to assign a volume to each particle, whose position within the box changes over time, and average that volume over the total period of the trajectory $\tau$. The way of partitioning space is to make the assignment that any point in space belongs to the particle which is closest to that point. For particle $i$, at any given time $t$:
$$ V_i (t) = \int_\text{box} \text{d}\textbf{r} \, \theta_i (\textbf{r}; t), \quad \text{where} \quad \theta_i (\textbf{r}; t) = \begin{cases} 1 \quad \text{if} \quad |\textbf{r}_i (t) - \textbf{r}| < |\textbf{r}_j (t) - \textbf{r}| \quad \forall \quad i \neq j, \\ 0 \quad \text{else}. \end{cases} $$
Here, $\textbf{r}_i (t)$ is the position of particle $i$ and $\textbf{r}_j (t)$ is the position of particle $j$. Then, the time average would be:
$$ \bar{V}_i = \frac{1}{\tau} \int_0^\tau \text{d} t \, V_i (t) $$
In practice, the time dependence is retrieved through finite sampling at equidistant time steps (I have $N_t$ lists of positions, each obtained a a given time).
I would like an efficient way to compute these values. Any good numerical approximation to the problem (I do not need extreme accuracy) would do. I don't know if the time part makes any difference, but perhaps someone can think of a way to make use of the time correlations.
At the moment my idea was to discretize the box into $N_p = N_x \times N_y \times N_z$ points and assign a discrete volume $\Delta V = \frac{V}{N_p}$ to each of them. Then loop over the list of particles to check which one is closest at each time and assign that $\Delta V$ to the closest particle. If the discretization is fine enough, this should be a good approximation to the integral, but I don't think it would be a very efficient method, especially given that I have many time steps to repeat this procedure.
Any ideas will be very welcomed.
In case anybody ends up reading this question, I found out that Voronoi tessellation does exactly what I want. I also found a C++ implementation online called Voro++ which I am successfully using for my calculations.