forked from DiedrichsenLab/surfAnalysis
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathsurf_resliceFS2WB.m
123 lines (106 loc) · 4.94 KB
/
surf_resliceFS2WB.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
function S=surf_resliceFS2WB(subj_name,subject_dir,atlas_dir,out_dir,varargin);
% functionsurf_resliceFS2WB(subj,subject_dir,atlas_dir,out_dir,varargin);
% Resampels a registered subject surface from freesurfer average to the new
% symmetric fs_LR_164 surface, standard in workbench.
% This allows things to happen exactly in atlas space - each vertex number
% corresponds exactly to a anatomical location
% For more information, see
% https://wiki.humanconnectome.org/download/attachments/63078513/Resampling-FreeSurfer-HCP_5_8.pdf
% INPUT:
% subj: subject name
% subjects_dir: freesurfer's SUBJECT_DIR
% atlas_dir: Directory with files for the standard_mesh (Atlas_templates/standard_mesh)
% out_dir: The new resampled files will be placed in (out_dir/subj)
% VARARGIN:
% 'hemisphere',[1 2] : left / right or both hemispheres
% 'align_surf': : Shift the surface to correct for freesurfer convention?
% 'surf_files',{'',''}: Surface files to be resampled
% {'.white','.pial','.inflated','.sphere.reg','.sphere'};
% 'curv_files',{'',''}: Curvature files to be resampled
% {'.curv','.sulc','.area'}
% ---------------------------
current_dir=pwd;
direct=[getenv('FREESURFER_HOME') filesep 'average' filesep 'surf' ];
if (isempty(subject_dir))
subject_dir=getenv('SUBJECTS_DIR');
end;
structname={'left','right'};
hem={'lh','rh'};
Hem={'L','R'};
hemisphere=[1:2]; % Do both hemispheres
surf_files={'.white','.pial','.inflated'};
curv_files={'.curv','.sulc','.area'};
smoothing=1;
align_surf=[1 1 1];
vararginoptions(varargin,{'smoothing','surf_files','curv_files','hemisphere','align_surf'});
% ----------------------------------------------------
% read freesurfer version
ver_file = fullfile(getenv('FREESURFER_HOME'),'build-stamp.txt');
if exist(ver_file)
fid = fopen(ver_file);
verStr = fgetl(fid);
fclose(fid);
verStr = strrep(verStr,'-v',' ');
verStr = textscan(verStr,'%s%f');
fsl_ver = verStr{2};
end
% ------------------------------------------------------
% create new output directory
new_dir = [out_dir filesep subj_name];
if (~exist(new_dir))
mkdir(new_dir);
end;
% -----------------------------------------------------
% Transform surface files to standard mesh
num_surf=length(surf_files);
cd(fullfile(subject_dir,subj_name,'surf'));
% -----------------------------------------------------
% Figure out the shifting of coordinate systems:
% Freesurfer uses vertex coordinates in respect to
% the center of the 256x256x256 image.
% Independent of the real zero point in the original image
% So to find a transform of the
% Mvox2surf: Transform of voxels in 256x256 image to surface vertices
% Mvox2space: Transform of voxel to subject space
anafile=fullfile(subject_dir,subj_name,'mri','brain.mgz');
[status,result] = system(['mri_info ' anafile ' --vox2ras-tkr']);
A=sscanf(result,'%f');
Mvox2surf=reshape(A,4,4)';
[status,result] = system(['mri_info ' anafile ' --vox2ras']);
A=sscanf(result,'%f');
Mvox2space=reshape(A,4,4)';
Msurf2space=Mvox2space*inv(Mvox2surf);
% -----------------------------------------------------
% Transform the surfaces from the two hemispheres
for h=hemisphere
% Convert reg-sphere
reg_sphere = [hem{h} '.sphere.reg.surf.gii'];
system(['mris_convert ' hem{h} '.sphere.reg ' reg_sphere]);
% -----------------------------------------------------
% Do all the surface files
for i=1:length(surf_files)
% Set up file names
file_name = [hem{h} surf_files{i} '.surf.gii'];
out_name = fullfile(new_dir,[subj_name '.' Hem{h} surf_files{i} '.164k.surf.gii']);
atlas_name = fullfile(atlas_dir,'resample_fsaverage',['fs_LR-deformed_to-fsaverage.' Hem{h} '.sphere.164k_fs_LR.surf.gii']);
% Convert surface to Gifti
system(['mris_convert ' hem{h} surf_files{i} ' ' file_name]);
system(['wb_command -surface-resample ' file_name ' ' reg_sphere ' ' atlas_name ' BARYCENTRIC ' out_name]);
A=gifti(out_name);
if (align_surf(i))
[A.vertices(:,1),A.vertices(:,2),A.vertices(:,3)]=surf_affine_transform(A.vertices(:,1),A.vertices(:,2),A.vertices(:,3),Msurf2space);
end;
save(A,out_name);
end;
% -----------------------------------------------------
% Do all the curvature files
for i=1:length(curv_files)
% Set up file names
file_name = [hem{h} curv_files{i} '.shape.gii'];
out_name = fullfile(new_dir,[subj_name '.' Hem{h} curv_files{i} '.164k.shape.gii']);
atlas_name = fullfile(atlas_dir,'resample_fsaverage',['fs_LR-deformed_to-fsaverage.' Hem{h} '.sphere.164k_fs_LR.surf.gii']);
% Convert surface to Gifti
system(['mris_convert -c ' [hem{h} curv_files{i}] ' ' [hem{h} surf_files{1}] ' ' file_name]);
system(['wb_command -metric-resample ' file_name ' ' reg_sphere ' ' atlas_name ' BARYCENTRIC ' out_name]);
end;
end;