Beauchamp:ALICE: Difference between revisions
No edit summary |
No edit summary |
||
Line 29: | Line 29: | ||
For subject YBK, the max intensity that worked best was 3500. Electrode volume was set to “2”, space was set to “1”. | For subject YBK, the max intensity that worked best was 3500. Electrode volume was set to “2”, space was set to “1”. | ||
The GUI then calls this script: | The GUI then calls this script: | ||
([tcsh -x 3dclustering.csh -CT_path ../CT/CT_highresRAI.nii -radius ' num2str(obj.settings.R) ' -interelectrode_space ' num2str(obj.settings.IS) ' -clip_value ' num2str(obj.settings.CV)]); | |||
This script contains the following commands: | This script contains the following commands: | ||
3dclust -savemask 3dclusters_r${r}_is${is}_thr${cv}.nii -overwrite -1Dformat -1clip $cv $is $r $ct > clst.1D | 3dclust -savemask 3dclusters_r${r}_is${is}_thr${cv}.nii -overwrite -1Dformat -1clip $cv $is $r $ct > clst.1D | ||
Line 44: | Line 44: | ||
echo "Clustered dataset saved as 3dclusters_r${r}_is${is}_thr${cv}.nii" | echo "Clustered dataset saved as 3dclusters_r${r}_is${is}_thr${cv}.nii" | ||
echo "Table of coordinates saved as clst.1D" | echo "Table of coordinates saved as clst.1D" | ||
Afni and SUMA were opened using this script: | |||
(['tcsh -x open_afni_suma.csh -CT_path ../CT/CT_highresRAI.nii -clust_set ' clust_set ' -clust_surf ' clust_surf ], '-echo') | |||
The script contains the following commands: | |||
The electrodes were manually selected using the 3-D clustering image in SUMA. Electrodes were selected individually and in order. SUMA was quit after selection was completed. | The electrodes were manually selected using the 3-D clustering image in SUMA. Electrodes were selected individually and in order. SUMA was quit after selection was completed. |
Revision as of 14:45, 21 July 2017
This page shows the steps used by the ALICE package to localize electrodes. These steps are called from within a Matlab GUI.
The GUI has three steps.
Here is the log file from the first step (an edited version of /Volumes/data/UT/YBK/ALICE/log_info/Step1_log.txt with comments added )
MRI scan selected: /Users/fosterlab/Documents/MATLAB/CTMR/DATA/YBK/ALICE/data/MRI/YBK_T1.nii
FS segmentation selected: /Users/fosterlab/Documents/MATLAB/CTMR/DATA/YBK/ALICE/data/FreeSurfer/t1_class.nii
This file was originally called YBK_fs_ribbon_rh_class.nii (renamed version of rh.ribbon.nii)
CT scan selected: /Users/fosterlab/Documents/MATLAB/CTMR/DATA/YBK/ALICE/data/CT/CT_highresRAI.nii
This file was originally called YBK_CT.nii
Aligning CT to MRI... This might take several minutes. Please wait... tcsh -x alignCTtoT1_shft_res.csh -CT_path CT_highresRAI.nii -T1_path../MRI/YBK_T1.nii
This shell script contains the following steps:
@Align_Centers -base $t1 -dset $ct 3dresample -input CT_highresRAI_shft.nii -prefix CT_highresRAI_res_shft.nii -master $t1 -dxyz 1 1 1 -rmode NN align_epi_anat.py -dset1 $t1 -dset2 CT_highresRAI_res_shft.nii -dset1_strip None -dset2_strip None -dset2to1 -suffix _al -feature_size 1 -overwrite -cost nmi -giant_move -rigid_body > status.txt 3dcopy CT_highresRAI_res_shft_al+orig CT_highresRAI_res_al.nii 3dcopy $t1 ./temp_ANAT.nii afni -com "SWITCH_UNDERLAY temp_ANAT.nii" -com "SWITCH_OVERLAY CT_highresRAI_res_al.nii"
The log file from the second step does not contain the called commands, the list below has been extracted from the Matlab file ctmrGUI.m
Three parameters are set in the GUI: electrode max intensity, electrode volume, inter electrode space.
For subject YBK, the max intensity that worked best was 3500. Electrode volume was set to “2”, space was set to “1”.
The GUI then calls this script:
([tcsh -x 3dclustering.csh -CT_path ../CT/CT_highresRAI.nii -radius ' num2str(obj.settings.R) ' -interelectrode_space ' num2str(obj.settings.IS) ' -clip_value ' num2str(obj.settings.CV)]);
This script contains the following commands:
3dclust -savemask 3dclusters_r${r}_is${is}_thr${cv}.nii -overwrite -1Dformat -1clip $cv $is $r $ct > clst.1D # make sure the clusters all show up in afni with distinct colors 3drefit -cmap INT_CMAP 3dclusters_r${r}_is${is}_thr${cv}.nii # now resample the clusters, erode, dilate and cluster again # this helps separate the clusters that overlap 3dresample -prefix temp_clusts_rs0.5 -overwrite -rmode NN -dxyz 0.5 0.5 0.5 -inset 3dclusters_r${r}_is${is}_thr${cv}.nii 3dmask_tool -dilate_inputs -1 +2 -prefix temp_clusts_rs0.5_de2 -overwrite -inputs temp_clusts_rs0.5+orig 3dclust -savemask 3dclusters_r${r}_is${is}_thr${cv}.nii -overwrite -1Dformat -1clip $cv $is $r temp_clusts_rs0.5_de2+orig > clst.1D 3drefit -cmap INT_CMAP 3dclusters_r${r}_is${is}_thr${cv}.nii rm temp_clusts*.HEAD temp_clusts*.BRIK* IsoSurface -isorois+dsets -mergerois+dset -autocrop -o_gii 3dclusters_r${r}_is${is}_thr${cv}.gii -input 3dclusters_r${r}_is${is}_thr${cv}.nii echo "Clustered dataset saved as 3dclusters_r${r}_is${is}_thr${cv}.nii" echo "Table of coordinates saved as clst.1D"
Afni and SUMA were opened using this script: (['tcsh -x open_afni_suma.csh -CT_path ../CT/CT_highresRAI.nii -clust_set ' clust_set ' -clust_surf ' clust_surf ], '-echo') The script contains the following commands:
The electrodes were manually selected using the 3-D clustering image in SUMA. Electrodes were selected individually and in order. SUMA was quit after selection was completed.
This script contains the following commands for electrode selection:
#get the xyz coordinate in the volume # really only need this for the afni_sphere case # could be a tiny bit faster without this check most of the time plugout_drive $NPB \ -com 'GET_DICOM_XYZ' \ -quit # have suma report its current surface label - which cluster DriveSuma $NPB -com "get_label" set clustindex = `tail -2 $sumasurf|head -1` #v2.0. old version=-1 $sumasurf` set xyzstr = `tail -1 $surfcoords` # output from plugout is of form "RAI xyz: x y z" # we can use just part of that # if we are using the exact location for a sphere, let's mark that here echo $electrode_i $clustindex $xyzstr[3-5] $afni_sphere >> $surfcoords_i set clustval = `echo $clustindex | sed 's/roi//' |sed 's/(I,T,B)R=//'|\ sed 's/(I,T,B)numeric=//'` # make 1/3 bright and add it to the list # this doesn't change the data in any way here # much smaller memory leak, 1000 iterations less that 256MB total for suma set roistr = `ccalc -form "roi%3.3d" $clustval` # commented lines only useful for checking if roi's have already been recolored # for coloring with white, gray or other constant color, skip these lines # use grep to be sure this cluster number is legit and to check status # don't put any commands beween the grep and status check # grep $roistr temproilist.txt # If first time cluster has been identified, recolor to 1/3 brightness # could make else condition, rebrighten electrode # if ($status) then set roirgb = `grep $roistr tempcmap.niml.cmap` set rgbout = ("1.0" "1.0" "1.0") #old version: `1deval -a "1D:$roirgb[1-3]" -expr a/3` # change niml file (redirect, then rename rather than in place with "sed -i ..." because of MacOS oddities) cat tempcmap.niml.cmap | sed "s/$roirgb/$rgbout 1 $roirgb[5] ${roistr}/" > tempcmap2.niml.cmap mv tempcmap2.niml.cmap tempcmap.niml.cmap # DriveSuma $NPB -com surf_cont -switch_dset temp_marked_clusters.1D.dset DriveSuma $NPB -com surf_cont -load_cmap tempcmap.niml.cmap # echo $roistr >> temproilist.txt # endif
The third and final step is electrode projection. The log file from the third step does not contain the called commands, the list below has been extracted from the Matlab file runMethod1.m. "Method 1" (Hermes et. al 2010), the subject name "YBK" and the hemisphere "right" were selected in the GUI. The grid settings were then manually created consisting of a label "G", the desired electrodes from the previous selection in step 2 "[1:32]", and grid dimensions "4,8".
The script calls the following commands: (Note that these commands are directly from the MATLAB script, as no .csh file exists for it)
%% NOTES; % This script uses freesurfer surface with electrode position extracted from % high res CT 3dclusters (using ctmr scirpt from Dora). %% 1.1) generate surface to project electrodes to hemisphere = obj.settings.Hemisphere; if strcmp(hemisphere, 'Left') hemi = 'l'; else hemi = 'r'; end if exist(['./results/' subject '_balloon_11_03.img'])==0 % if using freesurfer: get_mask_from_FreeSurfer(subject,... % subject name './data/FreeSurfer/t1_class.nii',... % freesurfer class file './results/',... % where you want to safe the file hemi,... % 'l' for left 'r' for right 11,0.3); % settings for smoothing and threshold %Visualize the surface with afni or mricron %saved as subject_balloon_11_03, where 11 and 0.3 are the chosen parameters. end %% 1.2) Convert DICOM to NIFTI + Coregistration + 3dClustering + electrode selection and sorting % coregister CT to anatomical MR using Afni and preserving CT resolution. % extract electrodes clusters using 3D clustering % extract CM using AFNI-SUMA plug-in. CM = importdata('./data/3Dclustering/electrode_CM.1D'); %remove repeated electrodes. [~, index] = unique(CM.data(:,4), 'last'); CM = CM.data(index,[1:4]); elecCoord = [-CM(:,1:2) CM(:,3)]; elecNum = CM(:,4); %check for empty rows and put NANs elecmatrix = zeros(elecNum(end),3); % create empty array auxk = 2; elecmatrix(1,:) = elecCoord(1,:); for k=2:length(elecNum) if elecNum(k)-elecNum(k-1)~= 1 elecmatrix(elecNum(k-1)+1:elecNum(k)-1,:) = nan; auxk = elecNum(k); elecmatrix(auxk,:) = elecCoord(k,:); auxk = auxk+1; else elecmatrix(auxk,:) = elecCoord(k,:); auxk = auxk +1; end end save([mypath 'CM_electrodes_sorted_all.mat'],'elecmatrix'); %% 4) Tansform coordinates to subject ANAT space: load([mypath 'CM_electrodes_sorted_all.mat']); Tmatrix = dlmread('./data/coregistration/CT_highresRAI_res_shft_al_mat.aff12.1D'); Tmatrix2 = dlmread('./data/coregistration/CT_highresRAI_shft.1D'); T = [reshape(Tmatrix',4,3) [0 0 0 1]']'; T2 = [reshape(Tmatrix2',4,3) [0 0 0 1]']'; T3 = T2*T; coord_al_anatSPM = []; coord_al_anat = []; for i = 1:size(elecmatrix,1) coord = [-elecmatrix(i,1:2) elecmatrix(i,3) 1]; coord_al_anat(i,:) = T3\coord' ; % = inv(T)*coord' coord_al_anatSPM(i,:) = [-coord_al_anat(i,1:2) coord_al_anat(i,3)]; end %%%%%%% SPM alignement: % T = [0.9979 0.0622 0.0154 -13.5301;... % -0.0629 0.9971 0.0428 -7.2991;... % -0.0127 -0.0437 0.9990 -19.1397;... % 0 0 0 1.0000]; % %coord_al_anatSPM = []; % %for i = 1:size(elecmatrix,1) % coord = [elecmatrix(i,1:3) 1]; %spm mode % coord_al_anatSPM(i,:) = inv(T)*coord' ; %end %%%%%%%%%%%%%%%%% % check result: % figure, % plot3(coord_al_anatSPM(:,1),coord_al_anatSPM(:,2),coord_al_anatSPM(:,3),'.','MarkerSize',20); hold on; % plot3(elecmatrix(:,1),elecmatrix(:,2),elecmatrix(:,3),'.r','MarkerSize',20); legend('aligned'); elecmatrix = coord_al_anatSPM(:,1:3); save([mypath 'CM_electrodes_sorted_all_aligned.mat'],'elecmatrix'); %% 5) project electrodes 2 surface system(['rm ./results/' subject '_singleGrid*']); load([mypath 'CM_electrodes_sorted_all_aligned.mat']); % electrodes2surf(subject,localnorm index,do not project electrodes closer than 3 mm to surface) % electrodes are projected per grid with electrodes2surf.m % in this example there were 7 grids % electrodes2surf( % 1: subject % 2: number of electrodes local norm for projection (0 = whole grid) % 3: 0 = project all electrodes, 1 = only project electrodes > 3 mm % from surface, 2 = only map to closest point (for strips) % 4: electrode numbers % 5: (optional) electrode matrix.mat (if not given, SPM_select popup) % 6: (optional) surface.img (if not given, SPM_select popup) % 7: (optional) mr.img for same image space with electrode % positions % automatically saves: % a matrix with projected electrode positions % a nifti image with projected electrodes % saved as electrodes_onsurface_filenumber_inputnr2 electrodes_path = [mypath 'CM_electrodes_sorted_all_aligned.mat']; surface_path = [mypath subject '_balloon_11_03.img']; anatomy_path = './data/FreeSurfer/t1_class.nii'; for g=1:size(obj.settings.Grids,2) grid = obj.settings.Grids{g}; %find comas: comas = strfind(grid,','); gridLabel = grid(1:comas(1)-1); gridEls = str2num(grid(comas(1)+1:comas(2)-1)); gridSize = str2num(grid(comas(2)+1:end)); if min(gridSize)==1 %strip parm(1) = 0; parm(2) = 2; elseif min(gridSize)==2 %small grid 2xN parm(1) = 4; parm(2) = 1; elseif min(gridSize)>2 %big grids: 3xN, 4xN, 6xN, 8xN parm(1) = 5; parm(2) = 1; else disp('! WARNING: Grid cannot have dimension 0. Please add grid again.'); %log str = get(obj.controls.txtLog, 'string'); if length(str)>=obj.settings.NUM_LINES str = str( (end - (obj.settings.NUM_LINES-1)) :end); end set(obj.controls.txtLog, 'string',{str{:},'>! WARNING: Grid cannot have dimension 0. Please add grid again.'}); loggingActions(obj.settings.currdir,3,' >! WARNING: Grid cannot have dimension 0. Please add grid again.'); end try [out_els,out_els_ind]=electrodes2surf_FreeSurfer(subject,parm(1),parm(2),gridEls,electrodes_path, surface_path, anatomy_path, mypath); catch parm(2) = 2; [out_els,out_els_ind]=electrodes2surf_FreeSurfer(subject,parm(1),parm(2),gridEls,electrodes_path, surface_path, anatomy_path, mypath); end elecmatrix(gridEls,:) = out_els; end %save as subject_singleGrid_projectedElectrodes_FreeSurfer_X_parm1_parm2, %where X is the number of the generated file (1 for the first file with %parameters X, Y and 2 for second file with same parameters), %and parm1 and parm2 the 2nd and 3rd of the projection function. %% 6) combine electrode files into one % save all projected electrode locaions in a .mat file save([mypath subject '_projectedElectrodes_FreeSurfer_3dclust.mat'],'elecmatrix'); % make a NIFTI image with all projected electrodes [output,els,els_ind,outputStruct]=... position2reslicedImage2(elecmatrix,anatomy_path); for filenummer=1:100 outputStruct.fname=[mypath subject '_projectedElectrodes_FreeSurfer_3dclust' int2str(filenummer) '.img' ]; if ~exist(outputStruct.fname,'file')>0 disp(['saving ' outputStruct.fname]); % save the data spm_write_vol(outputStruct,output); break end end %% 7) generate cortex to render images: hemisphere = obj.settings.Hemisphere; if strcmp(hemisphere, 'Left') % from freesurfer: in mrdata/.../freesurfer/mri gen_cortex_click_from_FreeSurfer(anatomy_path,[subject '_L'],0.5,[15 3],'l',mypath); display_view = [270 0]; % load cortex load([mypath subject '_L_cortex.mat']); save([mypath '/projected_electrodes_coord/' subject '_L_cortex.mat'],'cortex'); elseif strcmp(hemisphere, 'Right') % from freesurfer: in mrdata/.../freesurfer/mri gen_cortex_click_from_FreeSurfer(anatomy_path,[subject '_R'],0.5,[15 3],'r',mypath); display_view = [90 0]; % load cortex load([mypath subject '_R_cortex.mat']); save([mypath '/projected_electrodes_coord/' subject '_R_cortex.mat'],'cortex'); end %% 8) plot electrodes on surface % load electrodes on surface load([mypath subject '_projectedElectrodes_FreeSurfer_3dclust.mat']); % save final folder save([mypath '/projected_electrodes_coord/' subject '_projectedElectrodes_FreeSurfer_3dclust.mat'],'elecmatrix') ctmr_gauss_plot(cortex,[0 0 0],0); el_add(elecmatrix,'r',20); label_add(elecmatrix) loc_view(display_view(1), display_view(2));