ok so i decided to use setup_source_space instead, with ico5 and white as its arguments, and i called save_as_surface using that, and then converted the resulting .gii into glb, and it seemed to work, the vertex counts from .time.gii and the model were the same, so I can’t understand why the indices wouldn’t have correspondence
i didn’t quite understand the part about getting a source space mesh
my aim is to just reproduce the output of stc.plot as it is, with the optional given surface( white, inflated)
this is the code i’ve got to to make and write out the stc:
def setup_source_localization(spacing : str = "ico4") -> None:
global fsaverage_subjects_dir
global source_spaces
global bem_solution
global head_to_mri_trans
global forward_solution
fsaverage_subjects_dir = mne.datasets.fetch_fsaverage(
verbose=True
)
fsaverage_subjects_dir = str(fsaverage_subjects_dir)
# for surf in ['white', 'pial', 'inflated']:
# source_spaces[surf] = mne.setup_source_space(
# subject=fsaverage_subject,
# spacing=spacing,
# surface=surf,
# subjects_dir=fsaverage_subjects_dir,
# add_dist=False,
# verbose=True,
# n_jobs=-1
# )
print("-------------------------------------------------\nCreating source space")
source_spaces["white"] = mne.setup_source_space(
subject=fsaverage_subject,
spacing=spacing,
surface="white",
add_dist=False,
verbose=True,
n_jobs=-1
)
print("-------------------------------------------------\nCreated source space")
print("-------------------------------------------------\nCreating bem")
bem_model = mne.make_bem_model(
subject=fsaverage_subject,
ico=4,
conductivity=(0.3, 0.006, 0.3),
verbose=True
)
bem_solution = mne.make_bem_solution(
bem_model,
verbose=True
)
print("-------------------------------------------------\nCreated bem")
print("-------------------------------------------------\nCreating forward solution")
coreg = mne.coreg.Coregistration(
info=raw.info,
subject=fsaverage_subject,
fiducials="auto"
)
coreg.fit_fiducials(
verbose=True
)
head_to_mri_trans = coreg.trans
forward_solution = mne.make_forward_solution(
info=raw.info,
trans=head_to_mri_trans,
src=source_spaces["white"],
bem=bem_solution,
eeg=True,
meg=False,
mindist=5.0,
verbose=True
)
print("-------------------------------------------------\nCreated forward solution")
def make_stc(tmin : float = 0.0, tmax : float | None = None, SNR : float = 3.0, method : str = "dSPM") -> None:
global noise_covariance
global inverse_operator
global stc
lambda2 = 1.0 / (SNR * SNR)
if tmax is None:
tmax = raw.times[-1]
# ------------------------------------------------------------------
# Create noise covariance
# ------------------------------------------------------------------
# This is a pragmatic covariance estimate for continuous EEG.
#
# For a proper experimental source-localization pipeline, it is
# preferable to estimate the noise covariance from a baseline period
# or from empty-room / prestimulus data where appropriate.
events = mne.make_fixed_length_events(
raw,
duration=tmax-tmin
)
epochs = mne.Epochs(
raw,
events,
tmin=tmin,
tmax=tmax,
baseline=None,
preload=True,
verbose=False,
proj=True,
picks='all',
)
noise_covariance = mne.compute_covariance(
epochs, tmin=tmin, tmax=tmax, method=["empirical"], rank=None, verbose=True
)
evoked = epochs.average().pick("eeg")
del epochs
# ------------------------------------------------------------------
# Create inverse operator
# ------------------------------------------------------------------
inverse_operator = mne.minimum_norm.make_inverse_operator(
info=evoked.info,
forward=forward_solution,
noise_cov=noise_covariance,
loose=0.2,
depth=0.8,
)
# ------------------------------------------------------------------
# Apply inverse to continuous EEG
# ------------------------------------------------------------------
stc = mne.minimum_norm.apply_inverse(
evoked,
inverse_operator,
lambda2,
method=method,
pick_ori=None,
verbose=True,
)
def get_stc_dir_from_path(p : str) -> str:
return str(Path(p) / "source_estimates")
def get_stc_file_name() -> str:
return "stc-white"
def gifti_to_glb(input_file, output_file):
# Load GIFTI
img = nib.load(input_file)
vertices = None
faces = None
vertices = np.array([el.data for el in img.get_arrays_from_intent("NIFTI_INTENT_POINTSET")])
faces = np.array([el.data for el in img.get_arrays_from_intent("NIFTI_INTENT_TRIANGLE")])
vertices = vertices.reshape(-1, vertices.shape[-1])
faces = faces.reshape(-1, faces.shape[-1])
vertices = vertices.astype(np.float32)
faces = faces.astype(np.int64)
mesh = trimesh.Trimesh(
vertices=vertices,
faces=faces,
process=False
)
mesh.export(output_file, file_type="glb")
def write_stc_to_folder(p : str) -> None:
global stc
global source_spaces
from pathlib import Path
output_path = Path(get_stc_dir_from_path(p))
output_path.mkdir(parents=True, exist_ok=True)
# Already have white source space
file_name = str(output_path) + f"/{get_stc_file_name()}"
stc.save_as_surface(file_name, source_spaces["white"])
stc.save(file_name, overwrite=True)
gifti_to_glb(file_name + "-lh.gii", file_name + "-lh.glb")
gifti_to_glb(file_name + "-rh.gii", file_name + "-rh.glb")
and this is the code for reading a .time.gii file:
private byte[] GZipBase64ToUncompressedBytes(string gzipBase64)
{
// Step 1: Decode the Base64 text into the raw gzip-compressed bytes
byte[] gzipped = Convert.FromBase64String(gzipBase64);
// Step 2: Decompress with GZip to recover the original bytes
using var input = new MemoryStream(gzipped);
using var gzip = new ZLibStream(input, CompressionMode.Decompress);
using var output = new MemoryStream();
gzip.CopyTo(output);
return output.ToArray();
}
private float[][] LoadTimeCourseData(string path)
{
XDocument timeSeries = XDocument.Load(path);
var dataArrays = timeSeries.Root!.Elements("DataArray");
var datas = dataArrays.Select((el) =>
{
string compressedData = el.Element("Data").Value;
var uncompressed = GZipBase64ToUncompressedBytes(compressedData);
float[] data = new float[uncompressed.Length / 4];
for (int i = 0; i < uncompressed.Length; i += 4)
{
float val;
if (el.Attribute("Endian")!.Value == "LittleEndian")
{
val = BinaryPrimitives.ReadSingleLittleEndian(
uncompressed.AsSpan(i, 4));
}
else
{
val = BinaryPrimitives.ReadSingleBigEndian(
uncompressed.AsSpan(i, 4));
}
data[i / 4] = val;
}
return data;
}).ToArray();
Console.WriteLine($"{datas.Length}");
return datas;
}
