— Start with new version of Glasser atlas in fsaverage space. Downloaded from here: https://osf.io/rne5b/ filenames ?h.Glasser2016.mgz — convert to nii mri_convert ?h.Glasser2016.mgz ?h.Glasser2016.nii — get a surface reconstruction of the MNI Atlas 3dZeropad -prefix ./MNI_Atlas_zp.nii -RL 256 -AP 256 -IS 256 /Applications/AFNI/MNI152_T1_2009c+tlrc recon-all -all -sid MNI_Atlas -i MNI_Atlas_zp.nii @SUMA_Make_Spec_FS -NIFTI -sid MNI_Atlas — use Freesurfer to move from fsaverage to MNI mri_surf2surf —srcsubject fsaverage —trgsubject MNI_Atlas —sval atlasmgz/?h.Glasser2016.nii —tval MNI_Atlas/SUMA/?h.Glasser_HCP.gii —mapmethod nnf —hemi ?h cd MNI_Atlas/SUMA — use SUMA to resample the surface to std.141 SurfToSurf -i_gii std.141.?h.smoothwm.gii -i_gii ?h.smoothwm.gii -prefix res.141.lh -mapfile std.141.MNI_Atlas_?h.niml.M2M -dset ?h.Glasser_HCP.gii -output_params NearestNode mv res.141.?h.?h.Glasser_HCP.niml.dset res.141.?h.Glasser_HCP.niml.dset <- $hemi is in both prefix and dset from above, otherwise it will protest to overwriting the 1D file. Rename now so it doesn’t look silly, and so SUMA will load this later. — make a fatter ribbon so that all of the grey voxels get a label (I guess that’s a judgement call. The method for *not* doing that should be obvious) 3dcalc -a ?h.ribbon.nii -datum float -prefix ?h.fl.ribbon.nii -expr 'a' 3dmerge -1blur_fwhm 2 -doall -datum float -prefix ?h.blur.ribbon.nii ?h.fl.ribbon.nii — Use the new fat-ribbon @surf_to_vol_spackle -spec std.141.MNI_Atlas_?h.spec -surfA smoothwm -surfB pial -maskset ?h.blur.ribbon.nii’<0.25..2>’ -surfset res.141.?h.Glasser_HCP.niml.dset -mode -prefix ?h.Glasser_HCP.fill - Pad numbers so that left and right are separable 3dcalc -a lh.Glasser_HCP.fill.nii.gz -expr 'a' -prefix lh.MNI_Glasser_HCP.nii 3dcalc -a rh.Glasser_HCP.fill.nii.gz -expr 'a+(step(a)*1000)' -prefix rh.MNI_Glasser_HCP.nii - Combine into a single volume 3dcalc -a lh.MNI_Glasser_HCP+orig. -b rh.MNI_Glasser_HCP+orig. -expr 'max(a,b)' -prefix MNI_Glasser_HCP.nii - smooth 3dLocalstat -stat mode -nbhd 'SPHERE(2)' -prefix hcp_mode2 Glasser_HCP_labels.nii 3dcalc -a mode.MNI_Glasser_HCP.nii -b MNI_Glasser_HCP.nii \ -expr '(step(b)-step(a))*b + a' -prefix smooth.MNI_Glasser_HCP.nii \ -datum short - Add Atlas-y things 3drefit -space MNI smooth.MNI_Glasser_HCP.nii 3dcopy smooth.MNI_Glasser_HCP.nii MNI_Glasser_HCP_v1.0.nii.gz 3drefit -labeltable label_tables.txt MNI_Glasser_HCP_v1.0.nii.gz (open .HEAD file in a text editor) * atlas points were generated separately (see below) * copy paste Atlas points, adjusted to new numbers, fix character count * fixed a few stragglers * erased history ============= ============= *** making atlas locations *** ============= ============= (rh.HCP-MNI.annot is taken from here: https://figshare.com/articles/HCP-MMP1_0_projected_on_fsaverage/3498446) mris_divide_parcellation MNI_Atlas rh ./rh.HCP-MNI.annot 75 rh.75_div.annot mris_divide_parcellation MNI_Atlas lh ./lh.HCP-MNI.annot 75 lh.75_div.annot this gives us a new annotation with each of the original parcellations divided up into equal pieces about 75 mm^2 in area mri_annotation2label --subject MNI_Atlas --hemi ?h --ctab ./?h.75_ctab.txt --annotation 75_div export the annotation information into a text file (in matlab) [lh_temp.vol, lh_temp.M, lh_temp.mr_parms, lh_temp.volsz] = load_mgh('lh.w-g.pct.mgh'); [rh_temp.vol, rh_temp.M, rh_temp.mr_parms, rh_temp.volsz] = load_mgh('rh.w-g.pct.mgh'); % get template .mgz information from another dataset from this same subject [vertices lh_temp.vol ctab] = read_annotation('../label/lh.75_div.annot'); [vertices rh_temp.vol ctab] = read_annotation('../label/rh.75_div.annot'); % read in the annotation created above save_mgh(lh_temp.vol, 'lh.75_div.mgh', lh_temp.M, lh_temp.mr_parms); save_mgh(rh_temp.vol, 'rh.75_div.mgh', rh_temp.M, rh_temp.mr_parms); % save the annotation as a .mgh file !mri_convert lh.75_div.mgh lh.75_div.nii !mri_convert rh.75_div.mgh rh.75_div.nii % convert to .nii still_to_do = [1:180]; for n = 1:180 a=num2str(n); while length(a) < 3 a=['0' a]; end [~,roi_name] = system(['n=' a '; cat rh.75_ctab.txt | awk -v N=$n ''$1==N{print $2}'' ']); roi_name = roi_name(1:end-1); % for each ROI, find the right lines in the annotation output file [~,num_divs] = system(['n=' a '; cat rh.75_ctab.txt | grep ' roi_name ' | tail -n 1 | awk -F _div ''{print $2}'' | awk ''{print $1}'' ']); num_divs = str2num(num_divs); % figure out how many divisions this ROI was split into if num_divs/2 == round(num_divs/2) % even div1 = (num_divs/2); [~,div_label1] = system(['cat rh.75_ctab.txt | grep ' roi_name '_div' num2str(div1) ' | awk ''{print $1}'' ']); div_label1 = div_label1(1:end-1); div2 = (num_divs/2)+1; [~,div_label2] = system(['cat rh.75_ctab.txt | grep ' roi_name '_div' num2str(div2) ' | awk ''{print $1}'' ']); div_label2 = div_label2(1:end-1); if div_label1 + 1 ~= div_label2 continue end % find the middle two divisions, check that they make sense system(['mri_surfcluster --in rh.75_div.nii --subject MNI_Atlas2 --hemi rh --surf pial --centroid --thmin ' div_label1 ' --thmax ' div_label2 ' --sum surfclust_output_div/rh.p.' a '.txt --nofixmni']) system(['mri_surfcluster --in rh.75_div.nii --subject MNI_Atlas2 --hemi rh --surf white --centroid --thmin ' div_label1 ' --thmax ' div_label2 ' --sum surfclust_output_div/rh.w.' a '.txt --nofixmni']) % on both the pial surface and the white surface, find the location of the centroid voxel within the two central divisions still_to_do(n) = 0; % note that this ROI completed successfully else % odd div = (num_divs/2)+.5; [~,div_label] = system(['cat rh.75_ctab.txt | grep ' roi_name '_div' num2str(div) ' | awk ''{print $1}'' ']); div_label = div_label(1:end-1); system(['mri_surfcluster --in rh.75_div.nii --subject MNI_Atlas2 --hemi rh --surf pial --centroid --thmin ' div_label ' --thmax ' div_label ' --sum surfclust_output_div/rh.p.' a '.txt --nofixmni']) system(['mri_surfcluster --in rh.75_div.nii --subject MNI_Atlas2 --hemi rh --surf white --centroid --thmin ' div_label ' --thmax ' div_label ' --sum surfclust_output_div/rh.w.' a '.txt --nofixmni']) still_to_do(n) = 0; % same as above, but slightly simpler when there's an odd number of divisions end end for n = 1:180 a=num2str(n); while length(a) < 3 a=['0' a]; end % read in the output files created above, parsing for the centroid location [A,B] = system(['head -n 36 surfclust_output_div/lh.p.' a '.txt | tail -n 1 | awk ''{print $4}'' ']); disp(B) % I think this might be here to make sure there are the right number of lines in the file. [A,B] = system(['head -n 35 surfclust_output_div/lh.p.' a '.txt | tail -n 1 | awk ''{print $5}'' ']); [A,C] = system(['head -n 35 surfclust_output_div/lh.w.' a '.txt | tail -n 1 | awk ''{print $5}'' ']); Coords_lh(n,1) = -mean([str2num(B) str2num(C)]); [A,B] = system(['head -n 35 surfclust_output_div/lh.p.' a '.txt | tail -n 1 | awk ''{print $7}'' ']); [A,C] = system(['head -n 35 surfclust_output_div/lh.w.' a '.txt | tail -n 1 | awk ''{print $7}'' ']); Coords_lh(n,2) = mean([str2num(B) str2num(C)]); [A,B] = system(['head -n 35 surfclust_output_div/lh.p.' a '.txt | tail -n 1 | awk ''{print $6}'' ']); [A,C] = system(['head -n 35 surfclust_output_div/lh.w.' a '.txt | tail -n 1 | awk ''{print $6}'' ']); Coords_lh(n,3) = -mean([str2num(B) str2num(C)]); % find each coordinate as the average of the x, y, and z locations between calculating the centroid vertex % on the pial and the grey matter. It should (ideally) be the same vertex, so this should end up % giving us the mid-cortical location of that vertex. expr=['step(4-(x-' num2str(Coords_lh(n,1)) ')*(x-' num2str(Coords_lh(n,1)) ')-(y-' num2str(Coords_lh(n,2)) ')*(y-' num2str(Coords_lh(n,2)) ')-(z-' num2str(Coords_lh(n,3)) ')*(z-' num2str(Coords_lh(n,3)) '))']; pref = ['coord_sph/out.lh.' a '.nii']; p_list = findstr(expr, '--'); for p = length(p_list):-1:1 expr = [expr(1:p_list(p)-1) '+' expr(p_list(p)+2:end)]; end % build the pieces of the 3dcalc command to make a sphere centered around the coordinate location I found system(['3dcalc -a MNI_Glasser_HCP+tlrc -expr "' expr '" -prefix ' pref ' -overwrite']) % and execute that 3dcalc command clear expr pref p_list Coords_lh [A,B] = system(['head -n 35 surfclust_output_div/rh.p.' a '.txt | tail -n 1 | awk ''{print $5}'' ']); [A,C] = system(['head -n 35 surfclust_output_div/rh.w.' a '.txt | tail -n 1 | awk ''{print $5}'' ']); Coords_rh(n,1) = -mean([str2num(B) str2num(C)]); [A,B] = system(['head -n 35 surfclust_output_div/rh.p.' a '.txt | tail -n 1 | awk ''{print $7}'' ']); [A,C] = system(['head -n 35 surfclust_output_div/rh.w.' a '.txt | tail -n 1 | awk ''{print $7}'' ']); Coords_rh(n,2) = mean([str2num(B) str2num(C)]); [A,B] = system(['head -n 35 surfclust_output_div/rh.p.' a '.txt | tail -n 1 | awk ''{print $6}'' ']); [A,C] = system(['head -n 35 surfclust_output_div/rh.w.' a '.txt | tail -n 1 | awk ''{print $6}'' ']); Coords_rh(n,3) = -mean([str2num(B) str2num(C)]); expr=['step(4-(x-' num2str(Coords_rh(n,1)) ')*(x-' num2str(Coords_rh(n,1)) ')-(y-' num2str(Coords_rh(n,2)) ')*(y-' num2str(Coords_rh(n,2)) ')-(z-' num2str(Coords_rh(n,3)) ')*(z-' num2str(Coords_rh(n,3)) '))']; pref = ['coord_sph/out.rh.' a '.nii']; p_list = findstr(expr, '--'); for p = length(p_list):-1:1 expr = [expr(1:p_list(p)-1) '+' expr(p_list(p)+2:end)]; end system(['3dcalc -a MNI_Glasser_HCP+tlrc -expr "' expr '" -prefix ' pref ' -overwrite']); clear expr pref p_list Coords_lh % same thing for the other hemisphere end (back in bash) for n in `count -dig 3 1 180`; do for hemi in 'lh' 'rh'; do coord=`3dCM -mask MNI_Glasser_HCP+tlrc out.$hemi.$n.nii` echo $hemi $n $coord >> coord_list.txt done done (aaaaaand back to matlab) fid=fopen('atlas_points.txt', 'w'); fid_c = fopen('coord_sph/coord_list.txt', 'r'); %for n = 1:180 while 1 L = fgetl(fid_c) if L == -1 break end if strcmp(L(1:2), 'lh') hemi_add=0; H='L_'; else % rh hemi_add=1000; H='R_'; end L_nums = str2num(L(3:end)); n=L_nums(1); line = find(cell2mat(Names(:,1)) == n); % I think I copied and pasted these into the workspace by hand long_parcel_name = Names{line, 2}; parcel_name = Names{line, 2}; parcel_num = n+hemi_add; a=num2str(n); while length(a) < 3 a=['0' a]; end fprintf(fid, ['\n'], ... H, parcel_name, parcel_num, parcel_num, L_nums(2), L_nums(3), L_nums(4), H, long_parcel_name) end fclose(fid) fclose(fid_c) ============= ============= *** done making atlas locations *** ============= Yeah, I know. This was probably more work than it was worth, and probably not worth the wait. Oh, well. It's done now. ============= # make some pretty surfaces mkdir surfs_mode2 cd surfs_mode2 IsoSurface -isorois+dsets -o hcp_mode2.gii -input ../hcp_mode2_ribbon+tlrc. -Tsmooth 0.3 30 # show these surfaces echo "Use this command to see surfaces" echo "suma -onestate -i surfs_mode2/hcp_mode2*.gii &"