From d609c149c5a18efbdf6643b9428ce0b8fffa8fd2 Mon Sep 17 00:00:00 2001
From: "Peter J. Molfese" <pmolfese@gmail.com>
Date: Thu, 21 Dec 2023 11:21:55 -0500
Subject: [PATCH] [MRG][ENH]: Add Ability to export STC files as GIFTI (#12309)

Co-authored-by: Eric Larson <larson.eric.d@gmail.com>
Co-authored-by: Daniel McCloy <dan@mccloy.info>
---
 doc/changes/devel/12309.newfeature.rst |  1 +
 mne/source_estimate.py                 | 72 ++++++++++++++++++++++++++
 mne/tests/test_source_estimate.py      | 28 ++++++++++
 3 files changed, 101 insertions(+)
 create mode 100644 doc/changes/devel/12309.newfeature.rst

diff --git a/doc/changes/devel/12309.newfeature.rst b/doc/changes/devel/12309.newfeature.rst
new file mode 100644
index 00000000000..8e732044a8e
--- /dev/null
+++ b/doc/changes/devel/12309.newfeature.rst
@@ -0,0 +1 @@
+Add method :meth:`mne.SourceEstimate.save_as_surface` to allow saving GIFTI files from surface source estimates, by `Peter Molfese`_.
diff --git a/mne/source_estimate.py b/mne/source_estimate.py
index b2d197d7b2f..19b23da7d60 100644
--- a/mne/source_estimate.py
+++ b/mne/source_estimate.py
@@ -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
@@ -1584,6 +1585,77 @@ 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,
+            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
+        -----
+        .. versionadded:: 1.7
+        """
+        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()
+            darrays.append(
+                nib.gifti.gifti.GiftiDataArray(
+                    data=(s["rr"] * scale_rr).astype(np.float32),
+                    intent="NIFTI_INTENT_POINTSET",
+                    datatype="NIFTI_TYPE_FLOAT32",
+                )
+            )
+
+            # Make the topology DataArray
+            darrays.append(
+                nib.gifti.gifti.GiftiDataArray(
+                    data=s["tris"].astype(np.int32),
+                    intent="NIFTI_INTENT_TRIANGLE",
+                    datatype="NIFTI_TYPE_INT32",
+                )
+            )
+
+            # 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].astype(np.float32),
+                    intent="NIFTI_INTENT_POINTSET",
+                    datatype="NIFTI_TYPE_FLOAT32",
+                )
+                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.
 
diff --git a/mne/tests/test_source_estimate.py b/mne/tests/test_source_estimate.py
index be31fd1501b..ebe1a369e4d 100644
--- a/mne/tests/test_source_estimate.py
+++ b/mne/tests/test_source_estimate.py
@@ -248,6 +248,34 @@ def test_volume_stc(tmp_path):
             assert_array_almost_equal(stc.data, stc_new.data)
 
 
+@testing.requires_testing_data
+def test_save_stc_as_gifti(tmp_path):
+    """Save the stc as a GIFTI file and export."""
+    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."""