How to produce a 3D density map of a time-depenent system of particles? I have a time-dependent system of varying number of particles (~100k particles). In fact, each particle represents an interaction in a 3D space with a particular strength. Thus, each particle has $(X,Y,Z;w)$ which is the coordinate plus a weight factor between 0 and 1, showing the strength of interaction in that coordinate. Here, I have uploaded 10 real-time snapshots of the system, with particles are represented as reddish small dots; the redder the dot, the stronger the interaction is.
The question is: how one can produce a 3D (spatial) density map of these particles, preferably in Matlab or Origin Pro 9 or ImageJ? Is there a way to, say, take the average of these images based on the red-color intensity in ImageJ?
Since I have the numerical data for particles $(X,Y,Z;w)$ I can analyze those data in other software as well. So, you are welcome to suggest any other analytical approach/software.
 A: You're basically looking for a smoothing algorithm, something that takes a collection of points and turns it into a density. This can be done with the help of a kernel, a weighting function $K(u)$ that satisfies the following two conditions:
\begin{align}
\int_{-\infty}^{+\infty} K(u)\,du=1 \\
K(u)=K(-u)\quad\forall u
\end{align}
The Wikipedia article I link lists many common functions that satisfy these conditions. Be warned: the resulting distribution will depend on the kernel you choose. I wrote a Fortran 90 code that does this for 2D distributions for some work I did during the PhD; I haven't tried it, but it should be easily extensible to 3D.
From here, the basic idea is the following: You take your continuous space and discretize it (turn it into discrete blocks). Then you run through your list of $(x,y,z)$ tuples and find the closest grid-cell to your tuple. When you do that, loop around the 1st nearest neighbors and see of this tuple can 'leak' into the adjacent cells--the comparing distances (r < resolution below) effectively is your indicator function).
for each x,y,z in File
   iX <- closest(x)
   iY <- closest(y)
   iZ <- closest(z)
   for k=iZ-1,iZ+1
      for j=iY-1,iY+1
         for i=iX-1,iX+1
            r = sqrt((x-xgrid(i))^2 + (y-ygrid(j))^2 + (z-zgrid(k))^2)
            if r < resolution then
               map(i,j,k) = map(i,j,k) + kernel(r,resolution)
            end if
         end for
      end for
   end for
end for

where 


*

*closest is a function finds the closest cell to that particular x,y,z value

*xgrid (and others) are cell-centered vectors that hold the values of the discrete grid

*kernel is the aforementioned kernel function; note that u=r/resolution for comparing the code to the kernel page functions

*map is the 3D vector that stores the density at each cell (probably safe to initialize it to some low value, e.g., 0 or 1e-10)


Other thoughts I had on this:


*

*You'll need to account for your weighting in the above, which probably will just be  w*kernel, though I'm not entirely sure on that.  

*File output might also be tricky, as 3D density plots usually require some fancy formats; perhaps Matlab or Origin can handle such plots rather easily, I do not know.  

*Edge cases (where the data point is right at the boundary) will be slightly tricky with the iY-1 limits (some programs will wrap around to iYmax when negative numbers are reached), so you'll want to pick a domain that is completely outside the points.

