Bug fixes for running with SimNIBS 4

.gitignore

@@ -0,0 +1,2 @@
+
+.DS_Store

README.md

@@ -1,3 +1,33 @@
+# E-field Dosing
+## Dependencies and where to download them
+- **SimNIBS** and any additional dependencies (depending on version): https://simnibs.github.io/simnibs/build/html/installation/simnibs_installer.html
+  - The **SimNIBS example dataset** is found at https://simnibs.github.io/simnibs/build/html/dataset.html
+- **FreeSurfer**: https://surfer.nmr.mgh.harvard.edu/fswiki/DownloadAndInstall
+  - For Windows users, it is possible to download and use FreeSurfer through the Linux subsystem by following the instructions here: https://surfer.nmr.mgh.harvard.edu/fswiki/FS7_wsl
+## Checking the FreeSurfer integration with SimNIBS 4
+1. In the Linux subsystem / Terminal, navigate to the folder where the NifTI files are located in the SimNIBS example dataset
+2. In the Linux subsystem / Terminal, run `recon-all -subjid ernie -all -i org/ernie_T1.nii.gz -T2 org/ernie_T2.nii.gz`
+    1. This should generate a `fs_ernie` folder in the `SUBJECTS_DIR` you specified during the FreeSurfer Installation
+    2. If you get a permission denied error, add `export SUBJECTS_DIR=<path to acceessable folder>` to the same file you added `source $FREESURFER_HOME/SetUpFreeSurfer.sh` while setting up freesurfer
+3. In the Linux subsystem / Terminal, run `charm ernie org/ernie_T1.nii.gz org/ernie_T2.nii.gz --fs-dir <path to fs_ernie folder>`
+## Checking the TAP installation
+1. In MATLAB, Open “run_subj.m”
+3. In MATLAB, change the filepaths between `%%%%% CHANGE TAP PARAMETER HERE: Begin %%%%%` and `%%%%% CHANGE TAP PARAMETER HERE: End %%%%%` as needed
+    1. Set `Scale_Efield_to_Value` to the desired E-field value in V/m (can use 100 to test)
+    2. If not running on a Linux cluster, comment out the line starting with `setenv`
+4. Change `subjects/Ernie_charm` or the subjects subfolder corresponding to the appropriate pipeline to `subjects/ernie` and copy the `ernie` folder from the example dataset to `subjects/ernie`
+5. Copy `M1_mask_TMS_target.nii.gz` from `TAP/subjects/Ernie_headreco` into “subjects/ernie”
+6. In MATLAB, run `run_subj('ernie', 'subjects', 'M1', 'M1_mask_TMS_target.nii.gz', 1)`
+    1. This should generate a `ernie_M1` folder, a `ernie_M1_scirun` folder, and at least one `M1_HTX.X` folder
+## Quality control
+1. For the SimNIBS head mesh, see https://simnibs.github.io/simnibs/build/html/tutorial/advanced/fix_affine_registration.html for quality control
+2. Open `subjects/ernie/M1_HT0.0/ernie_TMS_optimize_MagVenture_Cool-B65.msh`
+    1. Check that the target is at the expected location
+    2. Check that the coil is oriented such that the primary E-field induced by the coil points into the target gyrus
+        - Many coils have arrows on them indicating the current direction
+        - For *monophasic* pulses, the primary E-field points in the *opposite* direction as the current
+        - For *biphasic* pulses, the primary E-field points in the *same* direction as the current
+
 # TAP: TMS_Targeting_And_Analysis_Pipeline
 Custom MATLAB code to determine an optimal scalp placement of a TMS coil on the computational model of the subject's head. (1) Optimal coil placement is chosen based on maximizing the TMS-induced electric field in a brain region of interest which is derived from individual fMRI peak activity and computed prospectively for neuronavigation-supported TMS experiments. (2) The code's functionality allows also for accuracy assessment of the coil placement recorded during and analyzed after a TMS session.  
 

matlab/get_target.m

@@ -198,7 +198,7 @@
         end
      end
     else
-      brain_voxels=GM_fromMesh;  
+      brain_voxels=find(GM_fromMesh);  
     end
 
     [gm_x gm_y gm_z] = ind2sub(size(GM_fromMesh),brain_voxels); 

matlab/importBS_mod.m

@@ -0,0 +1,188 @@
+% Modified version of David LK Murphy ([email protected])'s importBS("BStext")
+% bsData = importBS("BStext")
+% importBS is a function that converts brainsight text into a MATLAB 
+% structure. 
+% input: "BStext" is the filename or fullpath filename of the Brainsight
+% data text export.
+% output:"bsData" contians session information: Brainsight 
+% application info, samples, targets, and landmarks. 
+%
+% bsData Structure:
+% bsData.Info
+% bsData.Targets
+% bsData.Samples
+% bsData.PlannedLandmarks
+% bsData.SessionLandmarks
+% bsData.SessionName
+% 
+% All structure fields are in MATLAB table format, with the exception of 
+% "Info" which is a cell array. Variable names of each table correspond to
+% those found in Brainsight.
+%   
+% Sample data includes location, timestamp, and MEP values if recorded by
+% the Brainsight EMG system. The MEP/EMG part of the table in the "Samples"
+% field can be filled out with the recorded BrainVision EMG data.
+% Version: 0v2
+% Date: April 1, 2022
+% Author: David LK Murphy ([email protected])
+
+%% Setup working directory
+function bsData = importBS_mod(BStext)
+% BStext should be the full path to the file.
+rtrnPath=pwd;
+casa = pwd;
+if nargin<1||exist(BStext)==0
+    [BStext,BSpath] =uigetfile([casa filesep '*Brainsight_Export.txt']);
+    BStext = [BSpath BStext];
+end
+
+if isnumeric(BStext)||isempty(BStext)
+%     alertFig = uifigure;
+    warndlg('Import canceled: no file selected');
+return
+end
+firstLine = 8; % first line of session target data. Lines 1-7 are BS info.
+%% Set up the Import Options and import the data
+opts = delimitedTextImportOptions("NumVariables", 1);
+
+% Specify range and delimiter
+opts.DataLines = [1, inf];
+opts.Delimiter = "\t";
+
+% Specify column names and types
+% opts.VariableNames = ["Version12", "Var2", "Var3", "Var4", "Var5", "Var6", "Var7", "Var8", "Var9", "Var10", "Var11", "Var12", "Var13", "Var14", "Var15", "Var16", "Var17", "Var18", "Var19", "Var20", "Var21", "Var22", "Var23", "Var24", "Var25", "Var26", "Var27", "Var28", "Var29", "Var30", "Var31", "Var32", "Var33", "Var34", "Var35"];
+% opts.SelectedVariableNames = "Version12";
+% opts.VariableTypes = ["char", "char", "char", "char", "char", "char", "char", "char", "char", "char", "char", "char", "char", "char", "char", "char", "char", "char", "char", "char", "char", "char", "char", "char", "char", "char", "char", "char", "char", "char", "char", "char", "char", "char", "char"];
+
+% Specify file level properties
+% opts.ExtraColumnsRule = "ignore";
+% opts.EmptyLineRule = "read";
+
+bsData0 =readtable([BStext], opts);
+
+%% Convert to output type
+bsData0 = table2cell(bsData0);
+numIdx = cellfun(@(x) ~isnan(str2double(x)), bsData0);
+bsData0(numIdx) = cellfun(@(x) {str2double(x)}, bsData0(numIdx));
+bsData.Info = bsData0(1:7,1);
+targStart = find(contains({bsData0{:,1}}','Target Name'));
+sampStart = find(contains({bsData0{:,1}}','Sample Name'));
+plnStart = find(contains({bsData0{:,1}}','Planned Landmark Name'));
+slnStart = find(contains({bsData0{:,1}}','Session Landmark Name'));
+snmStart = find(contains({bsData0{:,1}}','Session Name'));
+%% Clear temporary variables
+clear opts bsData0
+
+%% Set up the Import Options and import the data
+opts = delimitedTextImportOptions("NumVariables", 13);
+
+% Specify range and delimiter
+opts.DataLines = [targStart+1 sampStart-1];
+opts.Delimiter = "\t";
+
+% Specify column names and types
+opts.VariableNames = ["TargetName", "LocX", "LocY", "LocZ", "m0n0", "m0n1", "m0n2", "m1n0", "m1n1", "m1n2", "m2n0", "m2n1", "m2n2"];
+opts.SelectedVariableNames = ["TargetName", "LocX", "LocY", "LocZ", "m0n0", "m0n1", "m0n2", "m1n0", "m1n1", "m1n2", "m2n0", "m2n1", "m2n2"];
+opts.VariableTypes = ["string", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double" ];
+
+% Specify file level properties
+opts.ExtraColumnsRule = "ignore";
+opts.EmptyLineRule = "read";
+
+% Specify variable properties
+opts = setvaropts(opts, ["TargetName"], "WhitespaceRule", "preserve");
+opts = setvaropts(opts, ["TargetName"], "EmptyFieldRule", "auto");
+opts = setvaropts(opts, ["LocX", "LocY","LocZ"], "TrimNonNumeric", true);
+opts = setvaropts(opts, ["LocX","LocY", "LocZ"], "ThousandsSeparator", ",");
+
+% Import the data
+bsData.Targets = readtable([BStext], opts);
+%% Clear temporary variables
+clear opts bsData0
+%% Set up the Import Options and import the data
+opts = delimitedTextImportOptions("NumVariables", 35);
+
+% Specify range and delimiter
+
+opts.DataLines = [sampStart+1, plnStart-1];
+opts.Delimiter = "\t";
+
+
+
+% Specify column names and types
+opts.VariableNames = ["SampleName", "SessionName", "Index", "AssocTarget", "LocX", "LocY", "LocZ", "m0n0", "m0n1", "m0n2", "m1n0", "m1n1", "m1n2", "m2n0", "m2n1", "m2n2", "DistToTarget", "TargetError", "AngularError", "TwistError", "StimPowerA", "StimPulseInterval", "StimPowerB", "Date", "Time", "CreationCause", "CrosshairsDriver", "Offset", "Comment", "EMGStart", "EMGEnd", "EMGRes", "EMGChannels", "EMGWindowStart", "EMGWindowEnd"];
+opts.VariableTypes = ["string", "string", "double", "string", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "datetime", "string", "string", "string", "double", "string", "double", "double", "double", "double", "double", "double"];
+
+% Specify file level properties
+opts.ExtraColumnsRule = "ignore";
+opts.EmptyLineRule = "read";
+
+% Specify variable properties
+opts = setvaropts(opts, ["SampleName", "SessionName", "AssocTarget", "Time", "CreationCause", "CrosshairsDriver", "Comment"], "WhitespaceRule", "preserve");
+opts = setvaropts(opts, ["SampleName", "SessionName", "AssocTarget", "Time", "CreationCause", "CrosshairsDriver", "Comment"], "EmptyFieldRule", "auto");
+opts = setvaropts(opts, "Date", "InputFormat", "yyyy-MM-dd");
+opts = setvaropts(opts, ["StimPowerA", "StimPulseInterval", "StimPowerB", "EMGStart", "EMGEnd", "EMGRes", "EMGChannels", "EMGWindowStart", "EMGWindowEnd"], "TrimNonNumeric", true);
+opts = setvaropts(opts, ["StimPowerA", "StimPulseInterval", "StimPowerB", "EMGStart", "EMGEnd", "EMGRes", "EMGChannels", "EMGWindowStart", "EMGWindowEnd"], "ThousandsSeparator", ",");
+
+% Import the data
+bsData.Samples  = readtable([BStext], opts);
+
+
+
+
+%% Clear temporary variables
+clear opts bsData0
+opts = delimitedTextImportOptions("NumVariables", 4);
+
+opts.DataLines = [plnStart+1 slnStart-1];
+opts.Delimiter = "\t";
+
+% Specify column names and types
+opts.VariableNames = ["PlannedLandmarkName", "LocX", "LocY", "LocZ"];
+    opts.SelectedVariableNames = ["PlannedLandmarkName", "LocX", "LocY", "LocZ"];
+opts.VariableTypes = ["string", "double", "double", "double" ];
+% Specify file level properties
+opts.ExtraColumnsRule = "ignore";
+opts.EmptyLineRule = "read";
+
+% Specify variable properties
+opts = setvaropts(opts, ["PlannedLandmarkName"], "WhitespaceRule", "preserve");
+opts = setvaropts(opts, ["PlannedLandmarkName"], "EmptyFieldRule", "auto");
+opts = setvaropts(opts, ["LocX","LocY", "LocZ"], "TrimNonNumeric", true);
+opts = setvaropts(opts, ["LocX","LocY", "LocZ"], "ThousandsSeparator", ",");
+
+% Import the data
+bsData.PlannedLandmarks = readtable([BStext], opts);
+
+%% Clear temporary variables
+clear opts
+
+%% Set up the Import Options and import the data
+opts = delimitedTextImportOptions("NumVariables", 6);
+
+% Specify range and delimiter
+opts.DataLines = [slnStart+1 snmStart-1];
+opts.Delimiter = "\t";
+
+% Specify column names and types
+opts.VariableNames = ["SessionLandmarkName", "SessionName", "Used", "LocX", "LocY", "LocZ"];
+opts.SelectedVariableNames = ["SessionLandmarkName", "SessionName", "Used", "LocX", "LocY", "LocZ"];
+opts.VariableTypes = ["string", "string", "string", "double", "double", "double"];
+
+% Specify file level properties
+opts.ExtraColumnsRule = "ignore";
+opts.EmptyLineRule = "read";
+
+% Specify variable properties
+opts = setvaropts(opts, ["SessionLandmarkName", "SessionName", "Used"], "WhitespaceRule", "preserve");
+opts = setvaropts(opts, ["SessionLandmarkName", "SessionName", "Used"], "EmptyFieldRule", "auto");
+opts = setvaropts(opts, "LocX", "TrimNonNumeric", true);
+opts = setvaropts(opts, "LocX", "ThousandsSeparator", ",");
+
+% Import the data
+bsData.SessionLandmarks = readtable([BStext], opts);
+
+%% Clear temporary variables
+clear opts
+cd(rtrnPath)
+end

matlab/make_brainsight_files.m

@@ -1,4 +1,4 @@
-function make_brainsight_files(subjects_folder, subject, target, A, outputfolder, ht, REVERSE_COIL_CURRENT, sep, coord_sys)
+function make_brainsight_files(subjects_folder, subject, target, A, outputfolder, ht, REVERSE_COIL_CURRENT, sep, coord_sys, varargin)
   point=A(1:3,4);
   if (~exist([subjects_folder sep subject sep outputfolder ]))
        mkdir([subjects_folder sep subject sep' outputfolder ]);
@@ -15,15 +15,20 @@ function make_brainsight_files(subjects_folder, subject, target, A, outputfolder
   fprintf(f,'# Units: millimetres, degrees, milliseconds, and microvolts\n');
   fprintf(f,'# Encoding: UTF-8\n');
   fprintf(f,'# Notes: Each column is delimited by a tab. Each value within a column is delimited by a semicolon.\n');
-  fprintf(f,'# Sample\tName\tSession\tName\tIndex\tAssoc. Target\tLoc. X\tLoc. Y\tLoc. Z\tm0n0\tm0n1\tm0n2\tm1n0\tm1n1\tm1n2\tm2n0\tm2n1\tm2n2\tDist. to Target\tTarget Error\tAngular Error\tTwist Error\tDate\tTime\tCreation Cause\tCrosshairs Driver\tOffset\tComment\n');
-  datum=char(datetime('now','Format','yyyy-MM-dd HH:mm:ss'));
-  datum_1=datum(1:strfind(datum,' ')-1);
-  datum_2=datum(strfind(datum,' ')+1:end);
   if (REVERSE_COIL_CURRENT==0)
-    fprintf(f,['Subject' '\t1\tSession\t2\t1\t' subject target 'HT' num2str(ht) 'mm' ' \t' num2str(point(1)) '\t'  num2str(point(2)) '\t' num2str(point(3)) '\t' num2str(A(1,1),4) '\t' num2str(A(2,1),4) '\t' num2str(A(3,1),4) '\t'  num2str(A(1,2),4) '\t' num2str(A(2,2),4) '\t' num2str(A(3,2),4) '\t' num2str(A(1,3),4) '\t' num2str(A(2,3),4) '\t' num2str(A(3,3),4) '\t' '0.0000\t0.0000\t0.0000\t0.0000\t' datum_1 '\t' datum_2 '.097\tButton\tMouse\t0.0000\t(null)\n']);
+      target_name = [subject target 'HT' num2str(ht) 'mm'];
   else
-    fprintf(f,['Subject' '\t1\tSession\t2\t1\t' subject  target 'HT' num2str(ht) 'mm' '_REVERSECOILCURRENT' ' \t' num2str(point(1)) '\t'  num2str(point(2)) '\t' num2str(point(3)) '\t' num2str(A(1,1),4) '\t' num2str(A(2,1),4) '\t' num2str(A(3,1),4) '\t'  num2str(A(1,2),4) '\t' num2str(A(2,2),4) '\t' num2str(A(3,2),4) '\t' num2str(A(1,3),4) '\t' num2str(A(2,3),4) '\t' num2str(A(3,3),4) '\t' '0.0000\t0.0000\t0.0000\t0.0000\t' datum_1 '\t' datum_2 '.097\tButton\tMouse\t0.0000\t(null)\n']);  
+      target_name = [subject target 'HT' num2str(ht) 'mm_REVERSECOILCURRENT'];
+  end
+  if ~isempty(varargin)
+      fprintf(f, '# Target Name\tLoc. X\tLoc. Y\tLoc. Z\n');
+      fprintf(f, [target_name sprintf('\t%.4f\t%.4f\t%.4f', varargin{1}) '\n']);
   end
+  fprintf(f,'# Sample Name\tSession Name\tIndex\tAssoc. Target\tLoc. X\tLoc. Y\tLoc. Z\tm0n0\tm0n1\tm0n2\tm1n0\tm1n1\tm1n2\tm2n0\tm2n1\tm2n2\tDist. to Target\tTarget Error\tAngular Error\tTwist Error\tDate\tTime\tCreation Cause\tCrosshairs Driver\tOffset\tComment\n');
+  datum=char(datetime('now','Format','yyyy-MM-dd HH:mm:ss'));
+  datum_1=datum(1:strfind(datum,' ')-1);
+  datum_2=datum(strfind(datum,' ')+1:end);
+  fprintf(f,['Subject 1\tSession 2\t1\t' target_name ' \t' num2str(point(1)) '\t'  num2str(point(2)) '\t' num2str(point(3)) '\t' num2str(A(1,1),4) '\t' num2str(A(2,1),4) '\t' num2str(A(3,1),4) '\t'  num2str(A(1,2),4) '\t' num2str(A(2,2),4) '\t' num2str(A(3,2),4) '\t' num2str(A(1,3),4) '\t' num2str(A(2,3),4) '\t' num2str(A(3,3),4) '\t' '0.0000\t0.0000\t0.0000\t0.0000\t' datum_1 '\t' datum_2 '.097\tButton\tMouse\t0.0000\t(null)\n']);
   fclose(f);
 
 disp(['wrote file: ' outputfolder sep subject '_' target '_hair' sprintf('%.1f',ht) 'mm.txt' ]);

optimization_check.m

@@ -0,0 +1,10 @@
+msh = mesh_load_gmsh4(fnamehead);
+mesh_show_surface(msh, 'showSurface', true)
+hold on
+l = scatter3(target(1), target(2), target(3), 'filled', 'DisplayName', 'Target Center');
+l = [l plot3([target(1) target(1)+5*target_direction(1)], ...
+    [target(2) target(2)+5*target_direction(2)], ...
+    [target(3) target(3)+5*target_direction(3)], ...
+    'LineWidth', 2, 'DisplayName', 'E-field Direction')];
+hold off
+legend(l)

retrospective_analysis.m

@@ -9,7 +9,7 @@
 subj='Ernie_headreco';
 subjects_folder='subjects';
 target_name='M1';
-outputfolder='M1_retrosim'
+outputfolder='M1_retrosim';
 
 BS_Trans_mat = read_brainsight_file([ subjects_folder sep subj sep target_name sep 'ernie_headreco_no-conform_M1_hair0.5mm.txt' ]);