Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[MRG][ENH]: Add Ability to export STC files as GIFTI #12309

Merged
merged 30 commits into from
Dec 21, 2023
Merged
Show file tree
Hide file tree
Changes from 24 commits
Commits
Show all changes
30 commits
Select commit Hold shift + click to select a range
b4d4bf4
add functions to export STC as GIFTI
pmolfese Dec 18, 2023
0e4a0a8
added more attribution
pmolfese Dec 18, 2023
7c5b675
move gifti export functionality to source_estimate.py
pmolfese Dec 19, 2023
42bdda6
Merge branch 'main' into exportstc_gifti
pmolfese Dec 19, 2023
c5ab7aa
Merge branch 'main' of https://github.com/mne-tools/mne-python into e…
pmolfese Dec 19, 2023
dab0639
devel doc update
pmolfese Dec 19, 2023
a0e19ab
Merge branch 'exportstc_gifti' of https://github.com/pmolfese/mne-pyt…
pmolfese Dec 19, 2023
1be89f1
add test to write and check save_as_surface
pmolfese Dec 19, 2023
5ed43f3
Merge branch 'main' into exportstc_gifti
pmolfese Dec 19, 2023
0e8960e
Reusing existing get_decimated_surfaces from source_space, missed tha…
pmolfese Dec 19, 2023
0994c0e
Merge branch 'exportstc_gifti' of https://github.com/pmolfese/mne-pyt…
pmolfese Dec 19, 2023
1784e66
styling fixes
pmolfese Dec 19, 2023
7eaafcc
verbose breakage
pmolfese Dec 19, 2023
070f9bb
Update doc/changes/devel.rst
pmolfese Dec 20, 2023
8b8f766
Update mne/source_estimate.py
pmolfese Dec 20, 2023
8613719
Update mne/source_estimate.py
pmolfese Dec 20, 2023
f21209e
Update mne/tests/test_source_estimate.py
pmolfese Dec 20, 2023
e10f5d7
TST: Ping
larsoner Dec 20, 2023
d97ab1e
Merge remote-tracking branch 'upstream/main' into exportstc_gifti
larsoner Dec 20, 2023
1587da1
STY: pre-commit
larsoner Dec 20, 2023
f205937
STY: Caps and dots
larsoner Dec 20, 2023
bca2773
FIX: Name
larsoner Dec 20, 2023
dffd5a9
FIX: No nest
larsoner Dec 20, 2023
11d77df
MAINT: DRY code
larsoner Dec 20, 2023
e6b7b0a
Update mne/tests/test_source_estimate.py
pmolfese Dec 20, 2023
86a3429
Update doc/changes/devel/12309.newfeature.rst
pmolfese Dec 20, 2023
a0df09d
Update mne/source_estimate.py
pmolfese Dec 20, 2023
8a0b922
Update mne/source_estimate.py
pmolfese Dec 20, 2023
d78b121
Apply suggestions from code review
larsoner Dec 21, 2023
c3d0075
Update mne/source_estimate.py
larsoner Dec 21, 2023
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions doc/changes/devel/12309.newfeature.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
Add method :class:`mne.SourceEstimate.save_as_surface` to allow saving GIFTI files from surface source estimates, by `Peter Molfese`_.
pmolfese marked this conversation as resolved.
Show resolved Hide resolved
74 changes: 74 additions & 0 deletions mne/source_estimate.py
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,7 @@
_ensure_src_subject,
_get_morph_src_reordering,
_get_src_nn,
get_decimated_surfaces,
)
from .surface import _get_ico_surface, _project_onto_surface, mesh_edges, read_surface
from .transforms import _get_trans, apply_trans
Expand Down Expand Up @@ -1584,6 +1585,79 @@ def in_label(self, label):
)
return label_stc

def save_as_surface(self, fname, src, *, scale=1, scale_rr=1e3):
"""Save a surface source estimate (stc) as a GIFTI file.

Parameters
----------
fname : path-like
Filename basename to save files as.
Will write anatomical gifti plus time series gifti for both lh/rh,
pmolfese marked this conversation as resolved.
Show resolved Hide resolved
for example ``"basename"`` will write ``"basename.lh.gii"``,
``"basename.lh.time.gii"``, ``"basename.rh.gii"``, and
``"basename.rh.time.gii"``.
src : instance of SourceSpaces
The source space of the forward solution.
scale : float
Scale factor to apply to the data (functional) values.
scale_rr : float
Scale factor for the source vertex positions. The default (1e3) will
scale from meters to millimeters, which is more standard for GIFTI files.

Notes
-----
Creates gifti files for source solution and time courses of STC.

pmolfese marked this conversation as resolved.
Show resolved Hide resolved
.. versionadded:: 1.7
pmolfese marked this conversation as resolved.
Show resolved Hide resolved
"""
nib = _import_nibabel()
_check_option("src.kind", src.kind, ("surface", "mixed"))
ss = get_decimated_surfaces(src)
assert len(ss) == 2 # should be guaranteed by _check_option above

# Create lists to put DataArrays into
hemis = ("lh", "rh")
for s, hemi in zip(ss, hemis):
darrays = list()
larsoner marked this conversation as resolved.
Show resolved Hide resolved
darrays.append(
nib.gifti.gifti.GiftiDataArray(
data=s["rr"] * scale_rr,
larsoner marked this conversation as resolved.
Show resolved Hide resolved
intent="NIFTI_INTENT_POINTSET",
datatype="NIFTI_TYPE_FLOAT32",
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

do we need to actually cast the data to float32 here? I read in nibabel docs that they don't actually enforce that the datatype matches the descriptor, but since we can make it match, seems like we should? (I'm assuming internally for us it's float64 right?)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'll have to test it. I'm checking against the non-Python GIFTI world since that's my use case. We ran into issues where nilearn was writing float64 NIFTI files and that didn't cooperate well with the AFNI/FSL/SPM/Freesurfer programs. Perhaps we add it as an option to match?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For now let's just cast to and use float32. We can add float64 later if we really need it but I'd assume YAGNI until someone hits a use case where float32 isn't enough

)
)

# Make the topology DataArray
darrays.append(
nib.gifti.gifti.GiftiDataArray(
data=s["tris"],
larsoner marked this conversation as resolved.
Show resolved Hide resolved
intent="NIFTI_INTENT_TRIANGLE",
datatype="NIFTI_TYPE_INT32",
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

similar comment as above about casting (here to int32, not float32)

)
)

# Make the output GIFTI for anatomicals
topo_gi_hemi = nib.gifti.gifti.GiftiImage(darrays=darrays)

# actually save the file
nib.save(topo_gi_hemi, f"{fname}-{hemi}.gii")

# Make the Time Series data arrays
ts = []
data = getattr(self, f"{hemi}_data") * scale
ts = [
nib.gifti.gifti.GiftiDataArray(
data=data[:, idx],
larsoner marked this conversation as resolved.
Show resolved Hide resolved
intent="NIFTI_INTENT_POINTSET",
datatype="NIFTI_TYPE_FLOAT32",
larsoner marked this conversation as resolved.
Show resolved Hide resolved
)
for idx in range(data.shape[1])
]

# save the time series
ts_gi = nib.gifti.gifti.GiftiImage(darrays=ts)
nib.save(ts_gi, f"{fname}-{hemi}.time.gii")

def expand(self, vertices):
"""Expand SourceEstimate to include more vertices.

Expand Down
28 changes: 28 additions & 0 deletions mne/tests/test_source_estimate.py
Original file line number Diff line number Diff line change
Expand Up @@ -248,6 +248,34 @@ def test_volume_stc(tmp_path):
assert_array_almost_equal(stc.data, stc_new.data)


@testing.requires_testing_data
pmolfese marked this conversation as resolved.
Show resolved Hide resolved
def test_save_stc_as_gifti(tmp_path):
"""Save the stc as a gifti file and export."""
pmolfese marked this conversation as resolved.
Show resolved Hide resolved
nib = pytest.importorskip("nibabel")
surfpath_src = bem_path / "sample-oct-6-src.fif"
surfpath_stc = data_path / "MEG" / "sample" / "sample_audvis_trunc-meg"
src = read_source_spaces(surfpath_src) # need source space
stc = read_source_estimate(surfpath_stc) # need stc
assert isinstance(src, SourceSpaces)
assert isinstance(stc, SourceEstimate)

surf_fname = tmp_path / "stc_write"

stc.save_as_surface(surf_fname, src)

# did structural get written?
img_lh = nib.load(f"{surf_fname}-lh.gii")
img_rh = nib.load(f"{surf_fname}-rh.gii")
assert isinstance(img_lh, nib.gifti.gifti.GiftiImage)
assert isinstance(img_rh, nib.gifti.gifti.GiftiImage)

# did time series get written?
img_timelh = nib.load(f"{surf_fname}-lh.time.gii")
img_timerh = nib.load(f"{surf_fname}-rh.time.gii")
assert isinstance(img_timelh, nib.gifti.gifti.GiftiImage)
assert isinstance(img_timerh, nib.gifti.gifti.GiftiImage)


@testing.requires_testing_data
def test_stc_as_volume():
"""Test previous volume source estimate morph."""
Expand Down
Loading