Aggregate points and surfaces in 4D

[jo-mueller]

Hi Robert @haesleinhuepf ,

I am seeing some issues with using the timeslicer on 4D points/surface data in napari. For instance, using the label_to_surface() function from napari-process-points-and-surfaces throws an error:

ValueError: Input volume should be a 3D numpy array.

which comes from the marching_cubes function under the hood. Here is a small example script to reproduce the error:

import napari
import napari_process_points_and_surfaces as nppas
# Make a blurry sphere
s = 100
data = np.zeros((s, s, s), dtype=float)
x0 = 50
radius = 15

for x in range(s):
    for y in range(s):
        for z in range(s):
            if np.sqrt((x-x0)**2 + (y-x0)**2 + (z-x0)**2) < radius:
                data[x, y, z] = 1.0

viewer = make_napari_viewer()
viewer.add_image(image)

segmentation = image > filters.threshold_otsu(image)
viewer.add_labels(segmentation)

surf = nppas.label_to_surface(segmentation.astype(int))
viewer.add_surface(surf)

When introspecting the call to marching_cubes within the time_slicer function it is also evident that the image is somehow still a 4D image.

[haesleinhuepf]

Hi Johannes @jo-mueller ,

yes, the time slicer at the moment only supports ImageData and LabelsData (see documentation). A fix could be straightforward, e.g. here:
https://github.com/haesleinhuepf/napari-workflows/blob/main/src/napari_workflows/_workflow.py#L595

See also this related issue:

• Support layers and/or layer metadata in workflows napari-workflows#2

Best,
Robert

[jo-mueller]

Hm...how would you identify points/surface data by array shapes?

A points data layer would have to be identified by being a 2D array with the second dimension being 4:

elif len(value.shape) == 2 and value.shape[-1]  == 4:
    layer = _get_layer_from_data(viewer, value)
    new_value  = value[value[:, 0] == current_timepoint, 1:]
    # I'm not sure whether breaking points down to 2D would make sense :)
    arguments[i] = new_value
    layer.metadata[CURRENT_TIME_FRAME_DATA] = new_value

A surface layer would have to be identified by being a tuple of length 2 or 3 with the above constraints for points:

elif isinstance(value, tuple) * len(value) in [2, 3] * value[0].shape[-1] == 4 

Would it be an option to allow napari.layers as inputs as well? From this it would be more evident what type of object is passed. I could see type checking by object shapes/lengths becoming a bit of a hassle.

[jo-mueller]

To get back this issue - I think a solution may be to fuse the functionalities of convert_to_file_backed_timelapse with the TimelapseConverter from napari-stress.

This part of the function:

for f in range(max_time):
    print("Processing frame", f)

    # go to a specific time point
    _set_timepoint(viewer, f)

does essentially the same as frame_by_frame in napari-stress. Since the results of that loop is simply a list of timepoint-results, using the TimeLapseConverter for this would be straightforward.

[jo-mueller]

Hi @haesleinhuepf ,

I added a function aggregate(layer : LayerInput, viewer: napari.Viewer = None) -> Layer: here. It gives me some weird results I don't quite understand if I use it for instance with this workflow on the napari-stress page.

It looks like the workflow manager is confused about having a 4D object as an endpoint to a workflow that handles 3D data when the time slider is dragged.