The model
I treat each particle as a hard sphere, with gaussian blur in 3D. Because of the lower z-resolution of the confocal microscopy, the blur is anisotropic. For particle $i$, I can write down its contribution ($\mathrm{S}_i$) to the the image as,
\[\mathrm{S}_i = \mathrm{I}_i \cdot (\mathrm{P}_i \ast \mathrm{G})\]where $\mathrm{I}$ is the intensity of the spehere, $\mathrm{P}$ is a “Platonic” particle without blur, and $\mathrm{G}$ is the gaussian kernel. A good algorithm to simulate the Platonic particle can be found in the SI of this nice PRX paper about particle tracking. As a result, the model of a confocal image ($\mathrm{M}$) can be written as
\[\mathrm{M} = \sum_i{\mathrm{S}_i}.\]In a more operational way, I can simulat a confocal image, knowing a collection of particle locations, the radii of these particles and their intensitiy values, following the pseudo-code,
for i in range(n_particles):
position = positions[i]
intensity = intensities[i]
radius = radii[i]
P = generate_single_sphere(radius)
box = get_sphere_container(positoin)
M[box] += P * intensity
apply_kernel(M, G)
The Derivative of the Model
If we are comparing the model with the experimental data, $\mathrm{E}$, we can take the squared difference as a measure of cost ($\mathrm{C}$). A good model should then be the model that minimise the difference, written as,
\[\begin{aligned} \mathrm{C} &= \sum_x\sum_y\sum_z(\mathrm{M} - \mathrm{E})^2 \cdot \mathrm{W}\\ &= \underbrace{\sum_x\sum_y\sum_z}_{\text{sum over all voxels}}\mathrm{D}^2 \cdot \mathrm{W} \end{aligned}\]where $\mathrm{D}$ is the difference between the experimental image and simulation, and $\mathrm{W}$ is the weight for different voxels. I choose $\mathrm{W} = \mathbf{1}$ and treat all pixel equally. Another option would be setting $\mathrm{W} \propto \mathrm{E}$ to stress the importance of bright voxels.
Suppose we want to use the criteria to find the particle locations, using algorithms like gradient descent, then we will need the derivative of $\mathrm{C}$ w.r.t. particle locations. The derivative w.r.t. the $x$ component of particle $i$ can be written as
\[\begin{aligned} \frac{\partial \mathrm{C}}{\partial x_i} &= \sum_x\sum_y\sum_z \mathrm{W} \cdot \frac{\partial \mathrm{D}^2} {\partial x_i}\\[1ex] &= \sum_x\sum_y\sum_z 2 \mathrm{DW} \cdot \frac{\partial \mathrm{D}} {\partial x_i}\\[1ex] &= \sum_x\sum_y\sum_z 2 \mathrm{DW} \cdot \frac{\partial \left[ \sum_j{\mathrm{I}_j \cdot (\mathrm{P}_j \ast \mathrm{G}}) - \mathrm{E} \right]} {\partial x_i}. \end{aligned}\]We know the experimental image would not be affected by the location of particle $i$. Therefore, $\partial\mathrm{E}/\partial x_i = 0$, and the derivative is finally,
\[\begin{aligned} \frac{\partial \mathrm{C}}{\partial x_i} &= \sum_x\sum_y\sum_z 2 \mathrm{DW} \cdot \frac{ \partial \sum_j{[\mathrm{I}_j \cdot (\mathrm{P}_j \ast \mathrm{G}}) ]} {\partial x_i} \\[1em] &= \sum_x\sum_y\sum_z 2 \mathrm{DW} \cdot \sum_j \left( \frac{ \partial \left[ \mathrm{I}_j \cdot (\mathrm{P}_j \ast \mathrm{G}) \right] } { \partial x_i } \right) \\[1em] &= \sum_x\sum_y\sum_z 2 \mathrm{DW} \cdot \sum_j \left( \frac{\partial \mathrm{I}_j}{\partial x_i} (\mathrm{P}_j \ast \mathrm{G}) + \mathrm{I}_j \frac{\partial (\mathrm{P}_j \ast \mathrm{G})}{\partial x_i} \right) \end{aligned}\]Exploiting the property of convolution operator
\[\frac{d}{dx}(f \ast g) = f \ast \frac{dg}{dx},\]we can further rewrite the differentiation as
\[\begin{aligned} \frac{\partial \mathrm{C}}{\partial x_i} &= \sum_x\sum_y\sum_z 2 \mathrm{DW} \cdot \sum_j \left( \frac{\partial \mathrm{I}_j}{\partial x_i} (\mathrm{P}_j \ast \mathrm{G}) + \mathrm{I}_j \cdot \left( \mathrm{P}_j \ast\frac{\partial \mathrm{G}}{\partial x_i} \right) \right) \\[1em] \end{aligned}\]since the intensity ($\mathrm{I}_j$) and shape ($\mathrm{P}_j$) of particle $j$ would not be affected by the location of particle $i$, the derivative can be simplified to the following,
\[\begin{aligned} \frac{\partial \mathrm{C}}{\partial x_i} &= \sum_x\sum_y\sum_z 2 \mathrm{DW} \; \left( \frac{\partial \mathrm{I}_i}{\partial x_i} \cdot (\color{teal}{\mathrm{P}_i} \ast \mathrm{G}) + \mathrm{I}_i \cdot \left(\color{teal}{\mathrm{P}_i} \ast \frac{\partial \mathrm{G}}{\partial x_i} \right) \right) \end{aligned}\]This equation can be used to evalulate the Jacobian matrix for the least–square fitting. It will be very fast because we only have to conserder a small sub–image, where $\mathrm{P}_i$ and its blurred version is non–zero. Formally, we can write something like,
\[\frac{\partial \mathrm{C}}{\partial x_i} = \underbrace{\color{teal}{ \sum_{x^\prime}\sum_{y^\prime}\sum_{z^\prime} }}_{\text{where P}_i \ast \mathrm{G} > 0 } 2 \mathrm{DW} \; \left( \frac{\partial \mathrm{I}_i}{\partial x_i} \cdot (\color{teal}{\mathrm{P}_i} \ast \mathrm{G}) + \mathrm{I}_i \cdot \left(\color{teal}{\mathrm{P}_i} \ast \frac{\partial \mathrm{G}}{\partial x_i} \right) \right).\]And the new summation range is much smaller.
Numerical Implementation
I implemented the cost function in my tracking code, and use the least square optimisation to update the particle locations.