forked from francopestilli/life_scripts
-
Notifications
You must be signed in to change notification settings - Fork 0
/
s_fe_make_mt_ips_rois.m
173 lines (142 loc) · 5.47 KB
/
s_fe_make_mt_ips_rois.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
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
function s_fe_make_mt_ips_rois(hemisphere)
%
% This script makes the MT and IPS rois used to find the tract reported in
% Pestilli et al., LIFE paper.
%
% Retinotopy folder: /biac3/wandell7/data/Retinotopy/
%
% ROIs can be created by using these two function in sequence:
% (1) fs_annotationToLabelFiles(fs_subject,annotationFileName,hemisphere,labelsDir)
% (2) fs_labelFileToNiftiRoi(fs_subject,labelFileName,niftiRoiName,hemisphere,regMgzFile,smoothingKernel)
%
% Copyright by Franco Pestilli Stanford University 2014
% Get the base directory for the data
subjects = {...
'FP_96dirs_b2000_1p5iso', ...
'KW_96dirs_b2000_1p5iso', ...
'MP_96dirs_b2000_1p5iso', ...
'HT_96dirs_b2000_1p5iso', ...
'JW_96dirs_b2000_1p5iso', ...
'KK_96dirs_b2000_1p5iso' ...
};
if notDefined('hemisphere'), hemisphere = {'lh','rh'};end
if notDefined('annotationFileName')
annotationFileName = {'aparc','aparc.a2009s'};
end
labelFileName = {'superiorparietal','MT'};
fsSubjectsDir = getenv('SUBJECTS_DIR');
for iSbj = 2%1:length(subjects)
if ~isdir(fullfile(fsSubjectsDir,subjects{iSbj}))
fs_subject = matchSubject2FSSUBJ(subjects{iSbj});
else fs_subject = subjects{iSbj};
end
fsSubjectDir = fullfile(fsSubjectsDir,fs_subject);
% Create all the necessary label files
for ia = 1:length(annotationFileName)
fs_annotationToLabelFiles(fs_subject,annotationFileName{ia},[],fsSubjectsDir);
end
% Extract the label files for the ROI we want and make nifti ROIs
for ih = 1:length(hemisphere)
for il = 1:length(labelFileName)
% Build a name for the nifti ROI and create it.
labelRoiName = sprintf('%s.%s.label',hemisphere{ih},labelFileName{il});
niftiRoiName = labelRoiName;
niftiRoiName(niftiRoiName=='.') = '_';
niftiRoiFullPath = fullfile(fsSubjectDir,'label', niftiRoiName);
fs_labelFileToNiftiRoi(fs_subject,labelRoiName,niftiRoiFullPath,hemisphere{ih},[],[],fsSubjectsDir);
end
end
end
end % Main function
function FS_SUBJECT = matchSubject2FSSUBJ(subject)
switch subject
case {'FP_96dirs_b2000_1p5iso'}
FS_SUBJECT = 'pestilli_test';
case {'KW_96dirs_b2000_1p5iso'}
FS_SUBJECT = 'weiner';
case {'MP_96dirs_b2000_1p5iso'}
FS_SUBJECT = 'lmperry';
case {'HT_96dirs_b2000_1p5iso'}
FS_SUBJECT = 'takemura';
case {'JW_96dirs_b2000_1p5iso'}
FS_SUBJECT = 'winawer';
case {'KK_96dirs_b2000_1p5iso'}
FS_SUBJECT = 'knk';
otherwise
keyboard
end
end
%%%%%%%%%%%%%%%%%%%%%%%
function [fh,sh] = makeBrainMap(fe,t1,slice,axLims,figName,saveDir)
% Make a map of the RMSE WITH and WITHOUT the fascicle:
coords = feGet(fe,'roi coords') + 1;
xform = feGet(fe,'xform img 2 acpc');
% Cross-validate RMSE
rmse = feGetRep(fe, 'vox rmse');
img = feReplaceImageValues(nan(feGet(fe,'map size')),rmse,coords);
maxr = 50;
% Make anifti file from the rmse
ni = niftiCreate('data',mbaNormalize(img,[0,1]), ...
'qto_xyz',xform, ...
'fname','rmse', ...
'data_type',class(img));
% Open a figure
fh = mrvNewGraphWin(figName);
% Show the anatomy with the overlay
sh = mbaDisplayOverlay(t1, ni, slice, [], 'hot');
axis(axLims)
saveMap(fh,figName,saveDir,nanmean(img(:)),nanmedian(img(:)),nanstd(img(:)),maxr)
end
%---------------------------------%
function saveMap(fh,figName,saveDir,M,m,SD,maxfd)
% This helper function saves two figures for each map and eps with onlythe
% axis and a jpg with only the brain slice.
% The two can then be combined in illustrator.
%
% First we save only the slice as jpeg.
set(gca,'fontsize',16,'ztick',[-20 0 20 40], ...
'xtick',[-50 -25 0 25 50], ...
'tickdir','out','ticklength',[0.025 0])
axis off
saveFig(fh,fullfile(saveDir,'maps',figName),'tiff')
saveFig(fh,fullfile(saveDir,'maps',figName),'png')
% Then we save the slice with the axis as
% eps. This will only generate the axis
% that can be then combined in illustrator.
axis on
grid off
title(sprintf('mean %2.2f | median %2.2f | SD %2.2f', ...
M,m,SD),'fontsize',16)
zlabel('Z (mm)','fontsize',16)
xlabel('X (mm)','fontsize',16)
cmap = colormap(hot(255));
colorbar('ytick',linspace(0,1,5),'yticklabel', ...
{linspace(0,1,5)*50}, ...
'tickdir','out','ticklength',[0.025 0],'fontsize',16)
saveFig(fh,fullfile(saveDir,'maps',figName),1)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function saveFig(h,figName,type)
% MAke sure the folder to save the figure exists
[p,f,e] = fileparts(figName);
[success,message] = mkdir(p);
if ~isempty(message), disp(sprintf('%s.',message));end
% Find out which type of figure and geenerate the proper printing command.
switch type
case {0,'jpeg','jpg'}
printCommand = (sprintf('print(%s, ''-djpeg90'',''-r500'' , ''-noui'', ''-opengl'', ''%s'')', num2str(h),figName));
case {1,'eps'}
printCommand = (sprintf('print(%s, ''-cmyk'', ''-depsc2'',''-tiff'',''-r500'' , ''-noui'', ''%s'')', num2str(h),figName));
case 'png'
printCommand = (sprintf('print(%s, ''-dpng'',''-r500'', ''%s'')', num2str(h),figName));
case 'tiff'
printCommand = (sprintf('print(%s, ''-dtiff'',''-r500'', ''%s'')', num2str(h),figName));
case 'bmp'
printCommand = (sprintf('print(%s, ''-dbmp256'',''-r500'', ''%s'')', num2str(h),figName));
otherwise
keyboard
end
% do the printing here:
fprintf('[%s] saving figure... \n%s\n',mfilename,figName);
eval(printCommand);
end