reorganize and clean up surface based minimum_norm

osl_ephys/report/src_report.py

@@ -76,16 +76,19 @@ def gen_html_data(config, outdir, subject, reportdir, logger=None, extra_funcs=N
     data["coregister"] = data.pop("coregister", False)
     data["beamform"] = data.pop("beamform", False)
     data["beamform_and_parcellate"] = data.pop("beamform_and_parcellate", False)
+    data["minimum_norm"] = data.pop("minimum_norm", False)
+    data["minimum_norm_and_parcellate"] = data.pop("minimum_norm_and_parcellate", False)
     data["fix_sign_ambiguity"] = data.pop("fix_sign_ambiguity", False)
 
     # Save info
-    if data["beamform_and_parcellate"]:
+    if data["beamform_and_parcellate"] or data['minimum_norm_and_parcellate']:
         data["n_samples"] = data["n_samples"]
     if data["coregister"]:
         data["fid_err"] = data["fid_err"]
-    if data["beamform_and_parcellate"]:
+    if data["beamform_and_parcellate"] or data['minimum_norm_and_parcellate']:
         data["parcellation_file"] = data["parcellation_file"]
-        data["parcellation_filename"] = Path(data["parcellation_file"]).name
+        if data["beamform_and_parcellate"]:
+            data["parcellation_filename"] = Path(data["parcellation_file"]).name
     if data["fix_sign_ambiguity"]:
         data["template"] = data["template"]
         data["metrics"] = data["metrics"]
@@ -242,15 +245,21 @@ def gen_html_summary(reportdir, logsdir=None):
     data["coregister"] = subject_data[0]["coregister"]
     data["beamform"] = subject_data[0]["beamform"]
     data["beamform_and_parcellate"] = subject_data[0]["beamform_and_parcellate"]
+    data["minimum_norm"] = subject_data[0]["minimum_norm"]
+    data["minimum_norm_and_parcellate"] = subject_data[0]["minimum_norm_and_parcellate"]
     data["fix_sign_ambiguity"] = subject_data[0]["fix_sign_ambiguity"]
 
     if data["coregister"]:
         subjects = np.array([d["filename"] for d in subject_data])
 
         fid_err_table = pd.DataFrame()
         fid_err_table["Session ID"] = [subject_data[i]["fif_id"] for i in range(len(subject_data))]
-        for i_err, hdr in enumerate(["Nasion", "LPA", "RPA"]):
-            fid_err_table[hdr] = [np.round(subject_data[i]['fid_err'][i_err], decimals=2) if 'fid_err' in subject_data[i].keys() else None for i in range(len(subject_data))]
+        if len(subject_data[0]['fid_err'])==4:
+            for i_err, hdr in enumerate(["Nasion", "LPA", "RPA", "Med(HSP-MRI)"]):
+                fid_err_table[hdr] = [np.round(subject_data[i]['fid_err'][i_err], decimals=2) if 'fid_err' in subject_data[i].keys() else None for i in range(len(subject_data))]   
+        else:
+            for i_err, hdr in enumerate(["Nasion", "LPA", "RPA"]):
+                fid_err_table[hdr] = [np.round(subject_data[i]['fid_err'][i_err], decimals=2) if 'fid_err' in subject_data[i].keys() else None for i in range(len(subject_data))]
         fid_err_table.index += 1 # Start indexing from 1
         data['coreg_table'] = fid_err_table.to_html(classes="display", table_id="coreg_tbl")
 

osl_ephys/source_recon/freesurfer_utils.py

@@ -93,20 +93,7 @@ def get_freesurfer_files(subjects_dir, subject):
         "std_brain_mri": op.join(os.environ["FREESURFER_HOME"], "subjects", "fsaverage", "mri", "T1.mgz"),
         "completed": op.join(surfaces_dir, "completed.txt"),
     }
-    #     "mni2mri_flirt_xform_file": op.join(surfaces_dir, "tranforms", "talairach.xfm"),
-    #     "mni_mri_t_file": op.join(surfaces_dir, "mni_mri-trans.fif"),
-    #     "bet_outskin_mesh_vtk_file": op.join(surfaces_dir, "outskin_mesh.vtk"),  # BET output
-    #     "bet_inskull_mesh_vtk_file": op.join(surfaces_dir, "inskull_mesh.vtk"),  # BET output
-    #     "bet_outskull_mesh_vtk_file": op.join(surfaces_dir, "outskull_mesh.vtk"),  # BET output
-    #     "bet_outskin_mesh_file": op.join(surfaces_dir, "outskin_mesh.nii.gz"),
-    #     "bet_outskin_plus_nose_mesh_file": op.join(surfaces_dir, "outskin_plus_nose_mesh.nii.gz"),
-    #     "bet_inskull_mesh_file": op.join(surfaces_dir, "inskull_mesh.nii.gz"),
-    #     "bet_outskull_mesh_file": op.join(surfaces_dir, "outskull_mesh.nii.gz"),
-    #     "std_brain": op.join(os.environ["FSLDIR"], "data", "standard", "MNI152_T1_1mm_brain.nii.gz"),
-    #     "std_brain_bigfov": op.join(os.environ["FSLDIR"], "data", "standard", "MNI152_T1_1mm_BigFoV_facemask.nii.gz"),
-    #     "completed": op.join(surfaces_dir, "completed.txt"),
-    # }
-
+
     # Coregistration files
     coreg_dir = op.join(fs_dir, "mne_src")
     os.makedirs(coreg_dir, exist_ok=True)
@@ -116,46 +103,10 @@ def get_freesurfer_files(subjects_dir, subject):
         "coreg_trans": op.join(coreg_dir, "coreg-trans.fif"),
         "coreg_html": op.join(coreg_dir, "coreg.html"),
         "source_space": op.join(coreg_dir, "space-src.fif"),
-        # "stc": op.join(coreg_dir, "stc-data", f"{subject}-lh.stc"),
-        # "filters_plot_cov": op.join(coreg_dir, "noise_cov.png"),
-        # "filters_plot_svd": op.join(coreg_dir, "noise_svd.png"),
+        "inverse_solution": op.join(coreg_dir, "{0}-inv.fif"),
+        "source_estimate_raw": op.join(coreg_dir, "src-raw"), # followed by -lh/rh.stc
+        "source_estimate_epo": op.join(coreg_dir, "src-epo"),
     }
-
-
-        # "info_fif_file": op.join(coreg_dir, "info-raw.fif"),
-        # "smri_file": op.join(coreg_dir, "scaled_smri.nii.gz"),
-        # "head_scaledmri_t_file": op.join(coreg_dir, "head_scaledmri-trans.fif"),
-        # "head_mri_t_file": op.join(coreg_dir, "head_mri-trans.fif"),
-        # "ctf_head_mri_t_file": op.join(coreg_dir, "ctf_head_mri-trans.fif"),
-        # "mrivoxel_scaledmri_t_file": op.join(coreg_dir, "mrivoxel_scaledmri_t_file-trans.fif"),
-        # "mni_nasion_mni_file": op.join(coreg_dir, "mni_nasion.txt"),
-        # "mni_rpa_mni_file": op.join(coreg_dir, "mni_rpa.txt"),
-        # "mni_lpa_mni_file": op.join(coreg_dir, "mni_lpa.txt"),
-        # "smri_nasion_file": op.join(coreg_dir, "smri_nasion.txt"),
-        # "smri_rpa_file": op.join(coreg_dir, "smri_rpa.txt"),
-        # "smri_lpa_file": op.join(coreg_dir, "smri_lpa.txt"),
-        # "polhemus_nasion_file": op.join(coreg_dir, "polhemus_nasion.txt"),
-        # "polhemus_rpa_file": op.join(coreg_dir, "polhemus_rpa.txt"),
-        # "polhemus_lpa_file": op.join(coreg_dir, "polhemus_lpa.txt"),
-        # "polhemus_headshape_file": op.join(coreg_dir, "polhemus_headshape.txt"),
-        # # BET mesh output in native space
-        # "bet_outskin_mesh_vtk_file": op.join(coreg_dir, "scaled_outskin_mesh.vtk"),
-        # "bet_inskull_mesh_vtk_file": op.join(coreg_dir, "scaled_inskull_mesh.vtk"),
-        # "bet_outskull_mesh_vtk_file": op.join(coreg_dir, "scaled_outskull_mesh.vtk"),
-        # # Freesurfer mesh in native space
-        # # - these are the ones shown in coreg_display() if doing surf plot
-        # # - these are also used by MNE forward modelling
-        # "bet_outskin_surf_file": op.join(coreg_dir, "scaled_outskin_surf.surf"),
-        # "bet_inskull_surf_file": op.join(coreg_dir, "scaled_inskull_surf.surf"),
-        # "bet_outskull_surf_file": op.join(coreg_dir, "scaled_outskull_surf.surf"),
-        # "bet_outskin_plus_nose_surf_file": op.join(coreg_dir, "scaled_outskin_plus_nose_surf.surf"),
-        # # BET output surface mask as nii in native space
-        # "bet_outskin_mesh_file": op.join(coreg_dir, "scaled_outskin_mesh.nii.gz"),
-        # "bet_outskin_plus_nose_mesh_file": op.join(coreg_dir, "scaled_outskin_plus_nose_mesh.nii.gz"),
-        # "bet_inskull_mesh_file": op.join(coreg_dir, "scaled_inskull_mesh.nii.gz"),
-        # "bet_outskull_mesh_file": op.join(coreg_dir, "scaled_outskull_mesh.nii.gz"),
-        # "std_brain": op.join(os.environ["FSLDIR"], "data", "standard", "MNI152_T1_1mm_brain.nii.gz"),
-
 
     # All Freesurfer files files
     files = {"surf": surf_files, "coreg": coreg_files, "fwd_model": op.join(fs_dir, "model-fwd.fif")}

osl_ephys/source_recon/minimum_norm.py

@@ -29,17 +29,12 @@
 def minimum_norm(
     outdir,
     subject,
-    preproc_file,
-    epoch_file,
+    data,
     chantypes,
     method,
     rank,
-    morph=False,
-    lambda2=0.1,
-    depth=0.8,
-    loose='auto',
-    freq_range=None,
-    pick_ori="normal",
+    depth,
+    loose,
     ):
     """Run minimum norm source localization.
     
@@ -49,10 +44,8 @@ def minimum_norm(
         Output directory.
     subject : str
         Subject ID.
-    preproc_file : str  
-        Preprocessed file.
-    epoch_file : str
-        Epoch file.
+    data : mnep.io.Raw, mne.Epochs  
+        Preprocessed data.
     chantypes : list
         List of channel types to include.
     method : str
@@ -61,77 +54,93 @@ def minimum_norm(
         Rank of the data covariance matrix.
     morph : bool, str
         Morph method, e.g. fsaverage. Can be False.
-    lambda2 : float
-        Regularization parameter.
     depth : float
         Depth weighting.
     loose : float
         Loose parameter.
-    freq_range : list
-        Band pass filter applied before source estimation.
-    weight_norm : str
-        Weight normalization.
-    pick_ori : str
-        Orientation to pick.
     reg : float
         Regularization parameter.
     reportdir : str
         Report directory.
-        
-        """
-
-    if preproc_file is None:
-        preproc_file = epoch_file
+    """
 
     log_or_print("*** RUNNING MNE SOURCE LOCALIZATION ***")
 
     fwd_fname = freesurfer_utils.get_freesurfer_files(outdir, subject)['fwd_model']
     coreg_files = freesurfer_utils.get_coreg_filenames(outdir, subject)
 
-    if epoch_file is not None:
-        data = mne.read_epochs(epoch_file, preload=True)
-    else:
-        data = mne.io.read_raw(preproc_file, preload=True)
-
-    # Bandpass filter
-    if freq_range is not None:
-        logger.info(f"bandpass filtering: {freq_range[0]}-{freq_range[1]} Hz")
-        data = data.filter(
-            l_freq=freq_range[0],
-            h_freq=freq_range[1],
-            method="iir",
-            iir_params={"order": 5, "ftype": "butter"},
-        )
-
-
-    if isinstance(data, mne.io.Raw):
-        data_cov = mne.compute_raw_covariance(data, method="empirical", rank=rank)
-    else:
-        data_cov = mne.compute_covariance(data, method="empirical", rank=rank)
     noise_cov = calc_noise_cov(data, rank, chantypes)
 
     fwd = mne.read_forward_solution(fwd_fname)
+    log_or_print(f"*** Making {method} inverse solution ***")
     inverse_operator = mne.minimum_norm.make_inverse_operator(data.info, fwd, noise_cov, loose=loose, depth=depth)
-    del fwd
+
+    log_or_print(f"*** Saving {method} inverse operator ***")
+    mne.minimum_norm.write_inverse_operator(coreg_files['inverse_operator'].format(method), inverse_operator, overwrite=True)
+    return inverse_operator
+
+
+def apply_inverse_solution(
+    outdir,
+    subject,
+    data,
+    method,
+    lambda2,
+    pick_ori,
+    inverse_operator=None,
+    morph="fsaverage",
+    save=False,
+    ):
+    """ Apply previously computed minimum norm inverse solution.
+    
+    Parameters
+    ----------
+    outdir : str
+        Output directory.
+    subject : str
+        Subject ID.
+    data : mne.io.Raw, mne.Epochs
+        Raw or Epochs object.
+    inverse_operator : mne.minimum_norm.InverseOperator
+        Inverse operator.
+    method : str
+        Inverse method.
+    lambda2 : float
+        Regularization parameter.
+    pick_ori : str
+        Orientation to pick.
+    morph : bool, str
+        Morph method, e.g. fsaverage. Can be False.
+    save : bool
+        Save source estimate (default: False).
+    """
+
+    coreg_files = freesurfer_utils.get_coreg_filenames(outdir, subject)
+    if inverse_operator is None:
+        inverse_operator = mne.minimum_norm.read_inverse_operator(coreg_files['inverse_operator'].format(method))
 
     log_or_print(f"*** Applying {method} inverse solution ***")
 
-    if epoch_file is not None:      
+    if isinstance(data, mne.Epochs):
         stc = mne.minimum_norm.apply_inverse_epochs(data, inverse_operator, lambda2=lambda2, method=method, pick_ori=pick_ori)    
     else:
         stc = mne.minimum_norm.apply_inverse_raw(data, inverse_operator, lambda2=lambda2, method=method, pick_ori=pick_ori)
 
     if morph:
+        log_or_print(f"*** Morphing source estimate to {morph} ***")
         src_from = mne.read_source_spaces(coreg_files['source_space'])
         morph = morph_surface(outdir, subject, src_from, subject_to=morph)
         morph.save(op.join(outdir, subject, "mne_src", morph), overwrite=True)
         stc = morph.apply(stc)
-
-    if epoch_file is not None:
-        stc.save(op.join(outdir, subject, "mne_src", "src-epo"), overwrite=True)
-    else:
-        stc.save(op.join(outdir, subject, "mne_src", "src-raw"), overwrite=True)
-
+
+    if save:     
+        log_or_print(f"*** Saving Source estimate ***")   
+        if isinstance(data, mne.Epochs):
+            stc.save(op.join(outdir, subject, "mne_src", "src-epo"), overwrite=True)
+        else:
+            stc.save(op.join(outdir, subject, "mne_src", "src-raw"), overwrite=True)
+    return stc
+
 
 def calc_noise_cov(data, data_cov_rank, chantypes):
     """Calculate noise covariance.
@@ -175,7 +184,7 @@ def calc_noise_cov(data, data_cov_rank, chantypes):
 
     bads = [b for b in data.info["bads"] if b in data_cov.ch_names]
     noise_cov = mne.Covariance(
-        noise_cov_diag, data_cov.ch_names, bads, data.info["projs"], nfree=1e10
+        noise_cov_diag, data_cov.ch_names, bads, data.info["projs"], nfree=data.n_times
     )
     return noise_cov
 

osl_ephys/source_recon/parcellation/parcellation.py

@@ -24,13 +24,15 @@
 from osl_ephys.utils.logger import log_or_print
 from osl_ephys.source_recon import freesurfer_utils
 
-def load_parcellation(parcellation_file):
+def load_parcellation(parcellation_file, subject=None):
     """Load a parcellation file.
 
     Parameters
     ----------
     parcellation_file : str
         Path to parcellation file.
+    subject : str
+        Subject ID. Only needed for FreeSurfer parcellations.
 
     Returns
     -------
@@ -40,9 +42,12 @@ def load_parcellation(parcellation_file):
 
     # check if it's a freesurfer parcellation
     if 'SUBJECTS_DIR' in os.environ:
-        avail = mne.label._read_annot_cands(os.path.join(os.environ["SUBJECTS_DIR"], 'fsaverage', 'label'))
+        if subject is None:
+            subject = "fsaverage"
+
+        avail = mne.label._read_annot_cands(os.path.join(os.environ["SUBJECTS_DIR"], subject, 'label'))
         if parcellation_file in avail:
-            labels = mne.label.read_labels_from_annot('fsaverage', parcellation_file)
+            labels = mne.label.read_labels_from_annot(subject, parcellation_file)
             if parcellation_file == 'aparc' or parcellation_file == "oasis.chubs":
                 labels = [l for l in labels if "unknown" not in l.name]
             elif parcellation_file == 'aparc.a2009s':
@@ -55,7 +60,7 @@ def load_parcellation(parcellation_file):
                 labels = [l for l in labels if "LOBE" in l.name] 
             return labels
 
-    # otherwise, load the parcellation file
+    # otherwise, load the nifti parcellation file
     parcellation_file = find_file(parcellation_file)
     return nib.load(parcellation_file)
 
@@ -363,7 +368,7 @@ def _get_parcel_timeseries(voxel_timeseries, parcellation_asmatrix, method="spat
     return parcel_timeseries, voxel_weightings, voxel_assignments
 
 
-def surf_parcellate_timeseries(subject_dir, subject, preproc_file, epoch_file, method, parcellation_file):
+def surf_parcellate_timeseries(subject_dir, subject, stc, method, parc):
     """Save parcellated data as a fif file.
     
     Parameters
@@ -372,23 +377,17 @@ def surf_parcellate_timeseries(subject_dir, subject, preproc_file, epoch_file, m
         Path to subject directory.
     subject : str
         Subject ID.
-    preproc_file : str
-        Path to preprocessed file.
-    epoch_file : str or None
-        Path to epoch file.
+    stc : mne.SourceEstimate
+        Source estimate.
     method : str
         Parcellation method. Can be 'pca_flip', 'max', 'mean', 'mean_flip', 'auto'
+    parc : str
+        Parcellation name.
     """
     fs_files = freesurfer_utils.get_freesurfer_files(subject_dir, subject)
-    if epoch_file is not None:
-        src_file = op.join(subject_dir, subject, "mne_src", "src-epo")
-    else:
-        src_file = op.join(subject_dir, subject, "mne_src", "src-raw")
 
-    labels = mne.read_labels_from_annot(subjects_dir=subject_dir, subject=subject, parc=parcellation_file)
-    labels = [l for l in labels if "unknown" not in l.name]
+    labels = load_parcellation(parc, subject=subject)
 
-    stc = mne.read_source_estimate(src_file, subject=subject)
     src = mne.read_source_spaces(fs_files['coreg']['source_space'])
     parcel_data = mne.extract_label_time_course(stc, labels, src, mode=method)
     return parcel_data