diff --git a/.gitignore b/.gitignore index e489273..1895afe 100644 --- a/.gitignore +++ b/.gitignore @@ -37,9 +37,9 @@ dist #ignore some random folders test_data -old + data_old .mypy_cache .pytest_cache -.ruff_cache +.ruff_cache \ No newline at end of file diff --git a/old/__pycache__/raw_cleaner.cpython-312.pyc b/old/__pycache__/raw_cleaner.cpython-312.pyc deleted file mode 100644 index 222fdaa..0000000 Binary files a/old/__pycache__/raw_cleaner.cpython-312.pyc and /dev/null differ diff --git a/old/calc_headmodels.py b/old/calc_headmodels.py deleted file mode 100644 index 36a5ff6..0000000 --- a/old/calc_headmodels.py +++ /dev/null @@ -1,135 +0,0 @@ -#%% -from cluster_jobs.meta_job import Job -#%% -import numpy as np -import mne -from mne.coreg import Coregistration -from os.path import join -import os -import pandas as pd -import pickle -from mne_bids import BIDSPath, read_raw_bids - -#%% -class HeadModelJob(Job): - - job_data_folder = 'headmodels' - data_file_suffix = '-trans.fif' - - def run(self, - subject_id, - base_dir, - savefig=True): - - #debug - #subject_id = 'CC320621' - base_path = '/home/schmidtfa/experiments/brain_age/data/data_cam_can' - mri_path = join(base_path, 'freesurfer') - - out_folder = join(base_path, 'headmodels', subject_id) - trans = 'fsaverage' - - - bids_path = BIDSPath(root=base_dir, subject=subject_id, session='rest', task='rest', - extension='.fif') - raw = read_raw_bids(bids_path=bids_path) - info = mne.pick_info(raw.info, mne.pick_types(raw.info, meg=True)) - - #%% do the coregistration - coreg = Coregistration(info, trans, mri_path) - coreg.set_scale_mode('3-axis') - coreg.fit_fiducials(verbose=True) - coreg.fit_icp(n_iterations=6, nasion_weight=2, verbose=True) - coreg.omit_head_shape_points(distance=5 / 1000) - coreg.fit_icp(n_iterations=20, nasion_weight=10, verbose=True) - - dists = coreg.compute_dig_mri_distances() * 1e3 # in mm - print(f"Distance between HSP and MRI (mean/min/max):\n{np.mean(dists):.2f} mm "f"/ {np.min(dists):.2f} mm / " - f"{np.max(dists):.2f} mm") - - #%% create and save a scaled copy of mri subject fs average - new_source_identifier = subject_id + '_from_template' - - #Dont forget this sdcales bem, atlas and source space automatically too - mne.coreg.scale_mri("fsaverage", new_source_identifier, - scale=coreg.scale, - subjects_dir=mri_path, - annot=True, overwrite=True) - - if os.path.isdir(out_folder) == False: - os.makedirs(out_folder) - file = open(join(out_folder, subject_id) + "info.pickle", 'wb') - pickle.dump(info, file) - file.close() - - mne.write_trans(join(out_folder, subject_id) + '-trans.fif', coreg.trans, overwrite=True) - print('Coregistration done!') - - if savefig: # REDO IT LATER - - # PLACEHOLDER: plot 3d stuff when xvfb is available - import matplotlib.pyplot as plt # Assuming 'info' is your data structure containing sensor and digitization information - head_mri_t = mne.transforms._get_trans(coreg.trans, "head", "mri")[0] - coord_frame = "head" - to_cf_t = mne.transforms._get_transforms_to_coord_frame(info, head_mri_t, coord_frame=coord_frame) - - sensor_locs = np.array([ch['loc'][:3] for ch in info['chs'] if ch['ch_name'].startswith('MEG')]) - sensor_locs = mne.transforms.apply_trans(to_cf_t['meg'], sensor_locs) - - # Extract Digitization Points (excluding fiducials) - head_shape_points = np.array([point['r'] for point in info['dig'] if point['kind'] == 4]) - head_shape_points = mne.transforms.apply_trans(to_cf_t['head'], head_shape_points) - - # Create a 2x2 subplot layout - fig = plt.figure(figsize=(10, 10)) - fig.suptitle(f'MEG - DIG Coregistration of {subject_id}', fontsize=16) - - # Axial view - ax1 = fig.add_subplot(221) - ax1.scatter(sensor_locs[:, 0], sensor_locs[:, 1], s=20, c='r', label='Sensors') - ax1.scatter(head_shape_points[:, 0], head_shape_points[:, 1], s=10, c='b', label='Head Shape') - ax1.set_title('Axial View') - ax1.set_xlabel('Distance (m)') - ax1.set_ylabel('Distance (m)') - ax1.legend() - - # Coronal view - ax2 = fig.add_subplot(222) - ax2.scatter(sensor_locs[:, 0], sensor_locs[:, 2], s=20, c='r') - ax2.scatter(head_shape_points[:, 0], head_shape_points[:, 2], s=10, c='b') - ax2.set_title('Coronal View') - ax2.set_xlabel('Distance (m)') - ax2.set_ylabel('Distance (m)') - - # Sagittal view - ax3 = fig.add_subplot(223) - ax3.scatter(sensor_locs[:, 1], sensor_locs[:, 2], s=20, c='r') - ax3.scatter(head_shape_points[:, 1], head_shape_points[:, 2], s=10, c='b') - ax3.set_title('Sagittal View') - ax3.set_xlabel('Distance (m)') - ax3.set_ylabel('Distance (m)') - - # 3D plot - ax4 = fig.add_subplot(224, projection='3d') - ax4.scatter(sensor_locs[:, 0], sensor_locs[:, 1], sensor_locs[:, 2], s=20, c='r', label='Sensors') - ax4.plot(sensor_locs[:, 0], sensor_locs[:, 1], sensor_locs[:, 2], color='k', linewidth=0.5) # Connect sensors with lines - ax4.scatter(head_shape_points[:, 0], head_shape_points[:, 1], head_shape_points[:, 2], s=10, c='b', label='Head Shape') - ax4.set_title('3D View') - ax4.grid(False) # Remove grid - ax4.axis('off') # Remove axis - - plt.tight_layout() - # save - plt.savefig(join(out_folder, subject_id) + '_coreg.png', dpi=300) - # not sohw ... plt.show() - - -# %% UNCOMMENT FOR TESTING -if __name__ == '__main__': - - subject_id = '19670901igsr' - - job = HeadModelJob(subject_id=subject_id) - job.run_private() - -# %% \ No newline at end of file diff --git a/old/preproc_sbg_raw.py b/old/preproc_sbg_raw.py deleted file mode 100644 index cc170c7..0000000 --- a/old/preproc_sbg_raw.py +++ /dev/null @@ -1,67 +0,0 @@ - - -from obob_mne.raw import Raw as RawTemplate -import numpy as np -import mne -from old.raw_cleaner import raw_cleaner - - -class Raw(RawTemplate): - sinuhe_root = '/home/schmidtfa/experiments/ncc/data/' - study_acronym = 'aw_ncc' - file_glob_patterns = ['%s_block%02d.fif', - '%s_block%d.fif'] - - - def run_cleaner(subject_id, - maxfilter = True, - ica = True, - ica_threshold = 0.5, - notch = False, - downsample_f=None, - l_pass = 40, - h_pass = 0.1): - - - n_blocks = Raw.get_number_of_runs(subject_id) - - #get average head pos - block_pos_l = [] - for block in np.arange(1, n_blocks + 1): - raw = Raw(subject_id, block_nr=block, preload=False) - block_pos_l.append(raw.info["dev_head_t"]['trans'][:3, 3]) - - blocks_pos = np.array(block_pos_l) - all_distances = np.sqrt(blocks_pos[:,0]**2 + blocks_pos[:,1]**2 + blocks_pos[:,2]**2) - mean_distance = np.median(all_distances) - mean_d_idx = (np.abs(all_distances - mean_distance)).argmin() - #TODO: maybe do np.median(np.array(block_pos_l), axis=0) - - raw_all, first_samples = [], [] - for block in np.arange(1, n_blocks + 1): - raw_tmp = Raw(subject_id, block_nr=block, preload=True) - - - #% Salzburg specific: make sure that if channels are set as bio that they get added correctly - if 'BIO003' in raw.ch_names: - raw.set_channel_types({'BIO001': 'eog', - 'BIO002': 'eog', - 'BIO003': 'ecg',}) - - mne.rename_channels(raw.info, {'BIO001': 'EOG001', - 'BIO002': 'EOG002', - 'BIO003': 'ECG003',}) - - #TODO: Only if trans average or trans default - eog_list = ['MEG0121','MEG0311' ,'MEG1211' ,'MEG1411'] - eog_list.extend([name for name in ['EOG001', 'EOG002'] if name in raw.ch_names]) - - # Determine the destination parameter for maxwell_filter - if trans_option == 'standard': - destination = (0, 0, 0.05) - elif trans_option == 'average' and input_files: - destination = compute_average_position(input_files) - else: - destination = None - - raw = raw_cleaner() \ No newline at end of file diff --git a/old/raw_cleaner.py b/old/raw_cleaner.py deleted file mode 100644 index f42f083..0000000 --- a/old/raw_cleaner.py +++ /dev/null @@ -1,60 +0,0 @@ -import mne -import numpy as np -from datetime import datetime -from os import listdir -from os.path import join - -#from rm_train import rm_train_ica - -def raw_cleaner(raw, - l_pass=None, - h_pass=None, - #ica part - ica=False, - ic_eog=True, - ic_ecg=True, - ic_muscle=False, - ic_train=False, - ic_train_freq=16.6, - ic_threshold=0.5, - #maxwell part - mw=True, - mw_coord_frame='head', - mw_destination=None, - mw_calibration_file = '../template_files/maxwell/sss_cal_sbg.dat', - mw_cross_talk_file = '../template_files/maxwell/ct_sparse_sbg.fif', - mw_st_duration=None,): - - - #maxfilter - if mw: - from utils.maxwell_utils import run_maxwell - maxwell_kwargs = {'coord_frame': mw_coord_frame, - 'destination': mw_destination, - 'calibration_file': mw_calibration_file, - 'cross_talk_file': mw_cross_talk_file, - 'st_duration': mw_st_duration} - - # find bad channels first - raw = run_maxwell(raw, **maxwell_kwargs) - - if np.logical_or(l_pass != None, h_pass != None): - raw.filter(l_freq=h_pass, h_freq=l_pass) - - #ica - if ica: - from utils.ica_utils import run_ica - ica_kwargs = {'resample_freq': 200, - 'eog': ic_eog, - 'ecg': ic_ecg, - 'muscle':ic_muscle, - 'train':ic_train, - 'train_freq': ic_train_freq, - 'ica_corr_thresh': ic_threshold} - - raw, ics = run_ica(raw, **ica_kwargs) - - return raw, ics - - else: - return raw \ No newline at end of file