read_tiff does not work for both non-imagej tiffs and greyscale tiffs

[BenCaiello]

Description
The read_tiff function seems like it will error for any non-imagej formatted image written by tifffile -- and perhaps any non-imagej tiff whatsoever -- as tifffile's read function will find the "imagej_metadata" attribute (which when not saved in imagej format = None) and that creates an error in NAVis, since NAVis expects the imagej_metadata attribute to not exist at all.

Secondarily, when I tried saving the image as with imagej metadata, because it was greyscale, navis still errored because of no support for greyscale / 2D tiffs (only 4D ZCYX tifs).

Is there a simple / alternate way to load numpy array data or .tiff image data into NAVis?

To Reproduce

import navis 
import tifffile
import numpy
directory = "C:/some/place/on/your/disc"
array = np.zeros([500,400])    ## should be any 2-dimensional grey-scale array which fails 
array = array[np.newaxis, np.newaxis,  :,  :]     ## trying to mimic "TCYX" order by adding axes changes nothing

#### Non-imageJ format fails:
tifffile.imwrite(directory, array)
navis.read_tiff(directory)

#### ImageJ format for greyscale / 2D fails:
tifffile.imwrite(directory, array.astype("uint8"), imagej = True, metadata = {})
navis.read_tiff(directory)
...

The arrays I'm trying to save are from an already skeletonized cell segmentations (hence the grayscale / 2D). What I've been working with tend to be small (smaller then the np.zeros faux-data in the example) because they are not necessarily neurons nor from very high resolution images. I have tried adding dimensions manually before saving to make them faux-4D, but somewhere in the save / load process they revert to 2D.

This error should apply to all non-imagej format tiffs:

This error should apply to greyscale imagej format tiffs:

Expected behavior
read_tiff reads tiffs

Your system

• NAVis version: 1.7.0
• • tifffile: 2024.2.12
• Python: 3.9
• Operating System: windows

[schlegelp]

Hi & thanks for the report! I will have a look at the issues with the header and channel.

Re 2d tiffs: navis.read_swc will (primarily) try to turn the image into a 3-dimensional neuron which is why I hadn't considered somebody wanting to read flat images. Could you perhaps elaborate a bit on what you would want to do with your images once they are in NAVis? It might also be helpful if you could share an example image :)

Edit:

Part of the problem is that TIFF is such a variable/flexible format when it comes to color vs greyscale, stacks vs single images, etc.

I'd also note that tiffile seems to parse headers/metadata from various different tools (imagej, ome, lsm and so on). I'm not sure we can cover all of these in navis.read_tiff off the bat but I'm happy to keep adding if people signal interest.

[BenCaiello]

Hello! Maybe I'm using NAVis improperly (trying 2D images), but I was hoping mainly just to get a diverse range of morphology measurements (like number of branches, branch points, etc.) and be able to measure those in 2 dimensions for any kind of cell type, like microglia or astrocytes, or others. I've not used NAVis before for its main purpose (.swc files and neurobiology), so maybe I'm just doing it wrong or have the wrong expectations haha.

Scikit-Image and other packages I've seen in python seem like they have much more limited morphology measuring abilities, although I had used skimage to skeletonize cell masks as preparatory step -- this is also why the actual tifffile's I'd been using are those written by tifffile from numpy arrays, after preparatory steps in skimage on cell segmentation masks.

I do think the fix for the primary error I encountered might be as easy changing how read_tiff parses imagj vs. non-imagej format -- in the navis/io/tiff_io.py file there is this block of code:

The issue seems like it might be that all .tiff files read by the tifffile package have this imagej_metadata attribute (even nonimagej formatted images), so the else statement (header = { } ) is never triggered. Instead, when the image is non-imagej, the header = None and you get the item assignment error I was encountering in line 179.
Maybe tifffile's behavior was changed / updated compared to when read_tiff was first written?

In any case I think a simple

if header is None:
       header = {}

in the right place should fix it without risking backward compatibility errors.

[schlegelp]

I think I have already fixed the header and channel errors in the fix_tiff branch. See this commit. It also introduces a couple of sanity checks such as checking that the data is actually 3D.

I'm holding off merging it in case we can make more changes to help you out. Let's assume we change read_tiff such that it accepts 2d images (we could e.g. simply fake the third dimension by reshaping the array). That would allow you to create either a VoxelNeuron (i.e. the image) or Dotprops (a point cloud). What you need is a TreeNeuron (=skeleton) to get things like cable length, branch points, etc. NAVis currently has limited capabilities to skeletonize Dotprops (see navis.skeletonize) but nothing for VoxelNeurons.

Since your data is already segmented, it would be best to go from the image straight to the skeleton. Implementing this has been on my TODO list for a while. You're welcome to have a crack at it yourself - shouldn't be too hard and I could provide initial pointers. Alternatively: it would be super helpful if you could share some example data for me to play with?

[schlegelp]

Actually, since your data is segmented, you should be able to generate a mesh using marching cubes, read that into NAVis and then skeletonize it. Again, would be helpful to get an idea of what your data looks like.

[BenCaiello]

Thanks for the help! I'm asking to see if I can share one of the data / images -- I work in a shared resource laboratory, so what I've been working with is not actually my / my lab's data -- or I might try to find a public brain / neural image dataset on zenodo and process that as an example.

What I already have successfully done with segmentation masks is:

1. Skeletonize the segmentations with scikit-image. Although perhaps this is not good or not necessary, since this may be it is better done in NAVis?

2. Iterate over each mask and extract it into a save-able form (numpy arrays) for export as .tiffs. These files are quite small because they only contain single cells / 2-dimensions. This step I'd done after skeletonization, but it could be done before as well.

If that code might be useful, I could share it, it's a pretty short bit of code to do those steps.

Also, I had the thought that maybe a new read_array function would be useful that accepts a numpy array + some metadata (like resolution), so that a .tiff intermediate is not necessary. But that might be more trouble for y'all than it is worth it.

How would marching cubes work in this case? I am not familiar with that algorithm.

[schlegelp]

That all sounds good. Just let me know if you have any such data that you can share.

For what it's worth: if you already have an array, you can just initialize a VoxelNeuron with it:

# With a dense array
arr = np.zeros((10, 10, 10))
arr[5, 5, 5] = 1
n = navis.VoxelNeuron(arr)

# With a sparse COO array
coo = np.arange(12).reshape(4, 3)
values = np.ones(4)
arr = np.hstack((coo, values.reshape(-1, 1)))
n = navis.VoxelNeuron(arr)

That said, I would recommend trying out marching cubes to generate a mesh - see e.g. scikit-learn's marching-cube implementation. It just takes a 3D array and turns it into vertices and faces for a mesh.

[BenCaiello]

OK, Thanks! I've managed to do the array --> marching cubes --> MeshNeuron --> TreeNeuron, although to get the marching cubes algorithm to work I had to concatenate my 2D image to itself (so that the z-axis had >1 length). Not sure how that might affect the output, I'm hoping not much. Now I'll explore what I can do in NAVis with it to get morphology measurements!

I also could do the array --> VoxelNeuron conversion easily, but wasn't sure what to do with the voxel neuron once it was created. It seems like the other way is better like you'd said.

[schlegelp]

Cool! Happy to help if you run into a dead end. Also: if this proves useful, perhaps you'd like to consider writing a tutorial for the documentation?

[schlegelp]

FYI, I added a new tutorial today that you may find useful: https://navis-org.github.io/navis/generated/gallery/0_io/zzz_tutorial_io_05_skeletonize/

[BenCaiello]

Thanks, sorry for the late response! I could consider writing up what I've done so far. It is very different from the tutorial you shared because my 2D .tiffs / segmentations comparatively are EXTREMELY small & sometimes minimal structure, since both being only in 2D prevents very large cells from being segmented & because are both neuronal & non-neuronal cells. In fact usually, the more straggly / interestingly shaped cells I've seen in the datasets I've worked with are astrocyte / microglia (not neuron), as axons are not well-defined.

So far, I've been able to get most of the morphology metrics I had hoping for like number of branch points, slab points, etc., using this code:

to_navis = image[k.bbox[0]:k.bbox[2],k.bbox[1]:k.bbox[3]]         ### This is a numpy array of the cell of interest                                  
                                         # extracted from a larger image using standard numpy slicing 
                                        # where k is the skimage regionprops() output for a given segmented cell
skeleton = ski.morphology.skeletonize(to_navis).astype('int')
to_navis = np.pad(skeleton[:,:,np.newaxis], 1)        ### marching cubes requires 3D with no degenerate dimensions
mesh = ski.measure.marching_cubes(to_navis, spacing = (1.0, 1.0, 0.01))
verts, faces, _, _ = mesh 
verts[:,2] = 0                      ### marching cubes output has variations in z-axis even though there really is no z-axis --> set all z's to 0
meshy_neury = nv.MeshNeuron((verts, faces))
tree_neuree = meshy_neury.skeletonize() 

The only issue I've run into is that in the conversion from mesh --> tree neuron I've sometimes seen straight skeletons (with no branches, visually speaking) have a node labeled as a branchpoint. Some of that had been z-axis variations (marching cubes must be padded into 3D to work, and I guess that algorithm can produce variations in z even across a 2D skeleton), but I fixed that by manually setting all the Z-values of the vertices of the mesh to 0 after running marching cubes.

Even with fixing the z-variations though, sometimes it visually looked like the algorithm that determines slab / branch points created a false branchpoint by cutting a corner. As in this very basic example:

from this segmented cell mask:

I get this skeleton:

Which can be converted into this mesh:

and finally this tree neuron:

But somehow this occurs:

Is this supposed to happened? I image in the vastly more complex 3D neuron skeletons I've seen in the tutorial you shared this kind of issue would not be easily seen, nor likely a big deal since there are so many branches compared to my itty-bitty 2D segmentations, haha.

[BenCaiello]

Hello!
Just wanted to check in on the status of this issue. It looked like there had been a PR & merge potentially working on some of these questions, but it'd be good to know where things are at.

[schlegelp]

Hi Ben. I did some re-work of read_tiff and other i/o functions. This includes trying to make it more robust against unexpected/missing file headers. It'd be great if you could update navis to 1.9.1 and see if your file is now read correctly.