DOC: First pass at SurfaceImage BIAP

[effigies]

Haven't fully covered what's in the HackMD, but did expand a bit where we made some decisions that didn't get written down. Also tried to pin down some terminology, as I think it's easy to speak imprecisely here. Some of this may be conventions that we need to enforce by transforming data on input (as we do with volumetric images come with affines that project into something besides RAS+).

Everything is still negotiable, so please feel free to make drastic suggestions.

Will start on implementations next week unless I feel inspired sooner.

cc @larsoner @htwangtw @NicolasGensollen @matthew-brett

[codecov]

Codecov Report

Merging #1056 (e278981) into master (a93f4b3) will decrease coverage by 0.32%.
The diff coverage is 42.85%.

@@            Coverage Diff             @@
##           master    #1056      +/-   ##
==========================================
- Coverage   92.22%   91.89%   -0.33%     
==========================================
  Files         100      100              
  Lines       12212    12322     +110     
  Branches     2136     2410     +274     
==========================================
+ Hits        11262    11323      +61     
- Misses        622      675      +53     
+ Partials      328      324       -4     

Impacted Files	Coverage Δ
nibabel/surfaceimages.py	42.85% <42.85%> (ø)
nibabel/streamlines/array_sequence.py	97.97% <0.00%> (-0.66%)	⬇️
nibabel/viewers.py	95.49% <0.00%> (-0.65%)	⬇️
nibabel/data.py	86.99% <0.00%> (ø)
nibabel/affines.py	100.00% <0.00%> (ø)
nibabel/arrayproxy.py	100.00% <0.00%> (ø)
nibabel/nicom/dicomwrappers.py	90.90% <0.00%> (ø)
nibabel/_h5py_compat.py
nibabel/streamlines/tractogram.py	99.68% <0.00%> (+<0.01%)	⬆️
nibabel/analyze.py	98.52% <0.00%> (+<0.01%)	⬆️
... and 13 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update a93f4b3...e278981. Read the comment docs.

[satra]

i wouldn't completely call cifti-2 a surface format. there is no geometry information stored in the file itself (unless someone hacked it. it's a point cloud that happens to have been extracted from some combination of a surface geometry and a volume geometry.

also many of the cifti-2 types are parcel based.

[htwangtw]

Indeed CIFTI supports multimodal data not purely for surface. It kind of deserve its own category.
IMO with the popularity of HCP data and fmriprep supporting output in fsLR template, surface based data support for CIFTI is too common to ignore. A major road block for user is data I/O. Working on how CIFTI-2 image relates to geometry template is a starting point none the less.

[effigies]

That's what I meant by "Pure data array". Fair point that it accepts data that has no geometric basis (parcels). Made a note of that case.

[NicolasGensollen]

Thanks @effigies !
I just made a quick pass but it looks good to me so far.

[bthirion]

Thx. Just left tiny comments.

[effigies]

@danjgale I would be curious about your thoughts, given that you have the surfplot package. What would make working with surface data smoother?

[danjgale]

@effigies thanks for the tag. TBH I found working with surfaces really easy when putting together surfplot, and I'm trying to recall any issues I had. When I eventually get back to polishing it up (soon🤞) and come across any issues, I can let you know.

The biggest challenges I've had primarily concern manipulating CIFTIs, such as changing the label table, creating a new CIFTI file from a modified data array, etc. Here's an example manipulation I've done in a different package, though I'm not sure how applicable/relevant these challenges are to the current topic.

[larsoner]

Just some preliminary small stuff, looks good so far!

[larsoner]

Might be worth detailing the number of dimensions supported in each of these, unless they all support arbitrary dimensions. Some or all (or none, I'm not sure!) of these could be a single scalar per vertex, so you can't for example have the n-dimensionality mentioned in the Data array definition above, you're stuck with n=0 (per vertex)

[effigies]

I considered that n=0, since one axis is the vertex.

[larsoner]

I think as long as there is a way to expose per-vertex data ndarray, we can just say "use NumPy-compatible tools" for this

[larsoner]

Having worked with VTK / Mayavi / PyVista / VisPy / mplot3d each a bit, I think this is going to be more difficult than it seems. This to me seems better tackled by a separate package, otherwise maintenance will be difficult.

I see the text below about NiBabel not necessarily providing each operation, so maybe adding to the Proposal below explicitly what functionality is out of scope would be good?

[larsoner]

I would put this at the top of this section -- when reading the list first, I assumed that the goal was to support/implement all prominent use cases.

[pep8speaks]

Hello @effigies, Thank you for updating!

• In the file nibabel/surfaceimages.py:

Line 112:38: W291 trailing whitespace

To test for issues locally, pip install flake8 and then run flake8 nibabel.

Comment last updated at 2022-02-11 02:14:46 UTC

[larsoner]

Sorry about missing the meeting today. I didn't get to the FreeSurfer stuff in the last two weeks but it looks like you started (yay!). Let me know if you want me to try it or add to it or anything else

[effigies]

Hey @larsoner, I've spent the morning fiddling with a collection API, if you want to have a look. I did the FreeSurfer stuff just to ground myself before diving into CIFTI and CaretSpecFile (now workbench specs), so probably won't be spending much more (self-directed) time on it, so you should feel free to play with it and see what it still needs to be useful.

[neurolabusc]

@effigies you may want to consider handling the mz3 format, the page includes Python, Blender, JavaScript code as well as a full specification and examples datasets. It attempts to live at the Pareto frontier of file size and speed. It is both faster and smaller than GIfTI. The layout design maps directly onto modern GPU architectures (and hence WebGL, Metal and Vulkan APIs). I developed this primarily for the sample files distributed with Surfice, to reduce download times, disk usage and the time to run the demos / regression testing. It is intentionally simple, and does not have the flexibility of other formats. It does not attempt to do everything, but it can fill a niche.

[effigies]

@neurolabusc Thanks! Opened an issue to avoid this getting lost...

[larsoner]

Just some preliminary testing from me, but looks like a good start. I can see the direction from the code and it makes sense.

[larsoner]

This is maybe going to be harder than we think:

>>> rr, tris = nib.freesurfer.read_geometry('/mnt/bakraid/data/structurals/fsaverage/surf/lh.white')
>>> from mne.surface import mesh_dist
>>> adjacency = mesh_dist(tris, rr)
>>> t0 = time.time(); dist_matrix = dijkstra(adjacency, directed=False, limit=10.); print(time.time() - t0)
...
numpy.core._exceptions.MemoryError: Unable to allocate 200. GiB for an array with shape (163842, 163842) and data type float64

So we'll need to improve SciPy's dijkstra to support sparse outputs start...

If we don't use the distance-based method and instead use the "smudging" operation by pushing data to neighboring vertices using a number of steps, this process is avoided. I guess this is a bit of an actual implementation detail, but we should think about these two methods as options with different tradeoffs.

[effigies]

Yes, full pairwise distances is painfully naive; not sure if iterated smudging is equivalent to distance-based (up to accumulated floating point error) but seems reasonable.

Mostly as a bookmark to myself, it might be worth looking into https://github.com/gilShamai/Fast-Pairwise-Geodesic-Distance-Computation to see if they have any good ideas.

[larsoner]

So far:

>>> from nibabel.tests.test_surfaceimages import FreeSurferSubject
>>> subj = FreeSurferSubject.from_subject('fsaverage', '/mnt/bakraid/data/structurals')
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/larsoner/python/nibabel/nibabel/tests/test_surfaceimages.py", line 210, in from_subject
    return klass.from_directory(Path(subjects_dir) / subject_id)
AttributeError: type object 'FreeSurferSubject' has no attribute 'from_directory'

So some unit tests are clearly needed :)

And if I change the call to use from_spec, then there is no self referred to in that class method

[larsoner]

This seems to work well as a subclass in local testing

[jbteves]

@larsoner @effigies where should I put some toy code as an example of converting an STC to a GIFTI? Should that go here, or should that be something we put in MNE first, and then add to this API later?

[effigies]

Would a gist be sensible in the short term?

[jbteves]

Sure, here you go:

https://gist.github.com/jbteves/384bfdd3777190f59620c3d988a314fd

[larsoner]

Feel free to merge and address my comments later @effigies , already looks like a great start!

[larsoner]

Or in syntactic sugar mode via __getitem__:

dlpfc_img = img[dlpfc_idcs]

and under the hood you just do the operation above?

[larsoner]

Add an example of trivial HDF5-based round-trip?

Native HDF5 format support
--------------------------
To support saving without loss of information, we also support serialization
to HDF5:
...

[larsoner]

One thing to note here is that betas will have collapsed over the time axis of the data, so the bold.header might be incorrect now, or otherwise need adjustment.

[larsoner]

Suggested change

	weights = normalize(gaussian(distances, sigma))
	weights = normalize(gaussian(distances, sigma)) # sparse matrix

[larsoner]

As you mentioned, this will kill most people's machines. Maybe you need to do this for GIfTI, but you shouldn't need to for the native HDF5 format, or for most operations (including plotting) -- it can just be done on the fly.

WDYT about an upsampler property or way to specify how to go from the subsampled / parcel space to the full-resolution space? By default / None could mean "map directly to the specified parcel, leave all else zero", but you in principle should be able to have this also be a scipy.sparse CSR or CSC matrix of shape (n_full_space, n_parcel) that you can @ with your data to get to the full-res space i.e., high-res surface for surfaces, or T1 / 1mm voxel spacing (usually) for volumetric.