3D point cloud with changing color in time

Hi,

I’m working in neuroscience, on 3d calcium recordings, and I’m trying to use napari to visualize the analysis process. I want to show 2 layers, with my stack in one, and each neuron activity after segmentation in the other. This means plotting points in 3D (x,y,z), and changing their colors with time.

I’m able to use napari to show my 4D (x,y,z,t) stacks, but I haven’t found a good way to show neuron activity.

Right now I’m doing the following (I have simplified for clarity) :

nb_points = < the number of neurons >
nb_t_steps = < the number of time steps >
points_coords = < an array with the (x,y,z) coordinates of my neurons
                              shape = [nb_points,3] >

viewer = napari.Viewer()

t = np.arange(nb_t_steps)
t = np.broadcast_to(t, (nb_points, nb_t_steps))
t = t.flatten(order="F")

Points = np.broadcast_to(
    points_coords, (nb_t_steps, nb_points, 3)
)
Points = np.reshape(
    Points, (nb_t_steps * nb_points, 3), order="C"
)

Points = np.c_[Points * 100, t]

act = points_activity.flatten()
color = np.c_[0.5 * (1 - act), 0.5 * (1 + act), 0.5 * (1 - act)]

viewer.add_points(Points, face_color=color, opacity=0.7)

This does exactly what I want, but it requires way to much memory because of all the broadcasting. So for a large number of neurons like 50k or more, it’s not a good solution.
Reading the napari documentation and forums, I haven’t found a satisfactory solution.

Does anyone have an idea of how to do this cleanly ?

Thanks a lot !

EmeEmu

Indeed you shouldn’t be forced to copy the x,y,z locations of all the points for all time just to change the colormap. I can’t off the top of my head think if we have an elegant way to solve this for the points layer yet (we probably should, but things will get more complicated). Maybe @kevinyamauchi has ideas.

We ran into something similar for our surface layer, see this example napari/surface_timeseries.py at e05e2efe1937b353b61e621b6aed93281be8b24f · napari/napari · GitHub and the data property of the Surface Layer - napari.layers.Surface - napari

Here you’re able to have a fixed surface, but changing values for each of the coordinates in time and we will handle the rest so you don’t need to broadcast.

We could think about doing something similar for the Points Layer. Curious how that sounds to you/ what others reading this think?

Suuuper interesting post, thanks for pushing the envelope @EmeEmu! :blush:

I have to agree with you @sofroniewn, I don’t see a clean API for this immediately, which is troublesome tbh :sweat_smile:. I think it’s likely a common use case.

Ideally, we should be able to broadcast changing properties over static points. But general properties seem more complicated to implement than just the values array for surfaces. I don’t know how I would expect to pass these values to napari in a way that was both intuitive and unambiguous.

The first workaround I would try is to have a 3D points layer with a callback on slider events that changes the properties on the points when the slider is moved. I’ll try to write a small example later today to add to our examples gallery. (I noticed in a live demo earlier this morning, because this is always when one notices these things :joy:, than viewer.dims.events.axis is gone, so I need to look up what the right event is now — @sofroniewn?) The nice thing about this is that (50k x time) is large but 50k is super small so this should be really fast and responsive.

1 Like

Ok thank you both @sofroniewn and @jni !
I’m very happy that I didn’t miss any obvious solution.

I understand that this feature would be hard to implement, but it would be extremely useful. It is probably the only missing feature for us to move our entire workflow to napari.

@jni , Thank you very much for adding your workaround to the example gallery, looking forward to it !

This is something that could one day we part of a custom properties widget. Say we did something like add_points(data_xyz, properties={'values': points_timeseries}), we could have a properties widget which allowed you to use a slider to map columns of the points_timeseries to the colormap. Not easy, but imaginable!!

See: Add example tying slider change to point properties change by jni · Pull Request #2582 · napari/napari · GitHub

Super interesting discussion, everyone. Thanks for the post, @EmeEmu !

I think @jni 's solution is probably the best solution given the current API. Looking forward, I could imagine expanding the layer properties to accept either arrays of shape (len(data),) or (len(data), n) for for n columns tied to some/all of the axes. Perhaps something like xarray would be nice so we can explicitly label the dims. Before this, I think we need to make a layer properties model to be able to store all of the state variables and update functions in a tidy manner.

Thanks you so much this workaround !! :smiley:

I copied your example directly to my machine, and I replaced 20 by 4000 and 300 by 50000 to test whether it would work with my data. At first it seemed to work perfectly, but when I clicked on the play button for the dimension 0 (time) slider, my RAM and Swap filled up in about 20 seconds (I have 16GB of RAM and 16GB of Swap). I tested this in more detail, and it seems that set_pts_properties saves some data each time it is called.

So for small arrays, this solution works wonders, but with more data I have the same problem I had before … it uses to much memory.

Another (less important question), If I comment out the image layer to show only the point layer, then napari thinks that my data is only 3D, and it doesn’t show a time slider. Is there a way to force napari to show the time slider ?

Again thanks a lot

Hmmm, that’s terrible and unexpected! I wonder actually whether it is related to dask and caching. Can you try adding from napari.utils import resize_dask_cache and then resize_dask_cache(0)?

At any rate, I’m confident we can track down the memory leak and plug it, don’t worry. :blush:

Yeah that’s an interesting one! Off the top of my head, nothing that isn’t a hack. I would suggest adding an invisible 4D layer containing e.g. 2 points, at the extremes of the range of the input data. But yeah, super hacky. :grimacing:

Hi,
So resize_dask_cache(0) did solve the problem !! No more filled RAM ! Thanks a lot ! :heart_eyes:

For future reference, using resize_dask_cache(0) does slow every down. When playing the slider, I need to reduce the fps to 2 in order to see all the frames, where 7 is enough when I don’t use resize_dask_cache(0). Even moving the slider by hand is extremely slow.

If the napari dev team does come up with a better solution in the future I think it would be a great improvement to an already amazing software. (should I create an issue on GitHub ?)
I will never the less tag your workaround as the solution !

again thanks a lot @jni

Nice. That seems to be a bug in dask, because the cache is by default 1/10th of RAM. I’ll try to make a napari-free example and raise an issue. But anyway, we’ll probably turn off caching by default in napari — that has been more trouble than it’s worth I think.

Yeah, two things here: first, the slowness is (I think/hope) all driven by the dask loading. So you should be able to make the image invisible and get your speed back. Do you observe that?

The other thing is the experimental asynchronous mode, which we’re still working on. You can turn it on with NAPARI_ASYNC=1 python myscript.py or by adding import os then os.environ['NAPARI_ASYNC'] = '1' at the top of your script (before importing napari). That should make things nice and speedy.

I think it’s probably a good idea to make an issue specifically to discuss the “points + time” API, thanks! I’ll handle the performance issues myself.

In the meantime, I’m glad you’re pretty much on your way. :blush:

OH! One more thing!

The issue with the timepoints being slow (I was surprised because each timepoint is pretty small — 16MB) is the default chunk size:

In [1]: image4d = da.random.random((4000, 32, 256, 256))

In [2]: image4d.chunksize
Out[2]: (80, 32, 64, 64)

So, without the cache, at each time point you are actually loading 80 timepoints! :man_facepalming:. (With the cache, you load the 80 timepoints once and then you reuse them until you reach the next chunk of 80.)

If you instead use:

image4d = da.random.random(
    (4000, 32, 256, 256),
    chunks=(1, 32, 256, 256),
)

your image chunks will be aligned with your slicing and everything will be nice and speedy. Does that pan out?

Ok I’ve updated the PR to use the bigger dataset (4000 timepoints, 50000 points), using dask for both the images and the points. Looks pretty cool I must say, even if it’s just a big ball pit. :joy: (click through for a short movie preview)

Yes this is much much better thank you so much !! :smiley:

So I open the following issue on GitHub.