view_sm(ROIs, type); Loads FS-FAST hemodynamic estimates and computes a selectivity map, displaying it on the slices. ROIs is a structure that contains information such as funcPath, nconditions, etc. Pre-defined selectivity ratios (given by the "type" parameter): 1) (positive - negative)/(positive + negative) in raw scanner units 2) (cond1 - cond2)/(hemodynamic offset) 3) (positive - negative)/(positive + negative) in percent signal change 4) abs(positive - negative)/(positive) in percent signal change $Id: view_sm.m,v 1.5 2003/08/27 16:07:03 nknouf Exp $
This function calls:
This function is called by:
0001 function view_sm(varargin) 0002 % view_sm(ROIs, type); 0003 % 0004 % Loads FS-FAST hemodynamic estimates and computes a selectivity 0005 % map, displaying it on the slices. 0006 % 0007 % ROIs is a structure that contains information such as funcPath, 0008 % nconditions, etc. 0009 % 0010 % Pre-defined selectivity ratios (given by the "type" parameter): 0011 % 0012 % 1) (positive - negative)/(positive + negative) in raw scanner units 0013 % 2) (cond1 - cond2)/(hemodynamic offset) 0014 % 3) (positive - negative)/(positive + negative) in percent signal change 0015 % 4) abs(positive - negative)/(positive) in percent signal change 0016 % 0017 % $Id: view_sm.m,v 1.5 2003/08/27 16:07:03 nknouf Exp $ 0018 0019 %% setup callback functions 0020 set(gcf,'KeyPressFcn', 'view_sm(''kbd'');'); 0021 set(gcf,'WindowButtonDownFcn', 'view_sm(''wbd'');'); 0022 0023 %% setup default values 0024 maxminfactor = 1.01; 0025 0026 %% run through the following if we are given 5 arguments 0027 if (nargin == 3) 0028 sm.ROIs = varargin{1}; 0029 sm.type = varargin{2}; 0030 sm.display = varargin{3}; 0031 0032 %% load the base image 0033 fprintf('view_sm: loading base volume...\n'); 0034 sm.base = fmri_ldbvolume([sm.ROIs.funcPath sm.ROIs.base '/' sm.ROIs.funcStem]); 0035 sm.base = sm.base(:,:,:,1); 0036 0037 %% load the hemodynamic volume 0038 fprintf('view_sm: loading hemodynamic volume...\n'); 0039 hemo_dir = sprintf('%s%s/h', sm.ROIs.funcPath, sm.ROIs.analysis); 0040 sm.hvol = fmri_ldbvolume(hemo_dir); 0041 0042 sm.cond1 = squeeze(sm.hvol(:,:,:, (sm.ROIs.cond1*2) + 1)); 0043 sm.cond2 = squeeze(sm.hvol(:,:,:, (sm.ROIs.cond2*2) + 1)); 0044 sm.cond1 = sm.cond1 + 1000; 0045 sm.cond2 = sm.cond2 + 1000; 0046 0047 blah = sm.cond1; 0048 save cond1.mat blah; 0049 blah = sm.cond2; 0050 save cond2.mat blah; 0051 0052 clear sm.hvol; 0053 0054 % load in hemodynamic offset volume 0055 fprintf('view_sm: loading hemodynamic offset volume...\n'); 0056 hemo_offset_path = sprintf('%s%s/h-offset', sm.ROIs.funcPath, sm.ROIs.analysis); 0057 sm.hoffsetvol = fmri_ldbvolume(hemo_offset_path); 0058 0059 switch (sm.type) 0060 case 1, 0061 % compute selectivity map 0062 sm.sm = (sm.cond1 - sm.cond2)./(sm.cond1 + sm.cond2); 0063 case 2, 0064 sm.sm = (sm.cond1 - sm.cond2)./(sm.hoffsetvol); 0065 case 3, 0066 sm.cond1 = sm.cond1 - 1000; 0067 sm.cond2 = sm.cond2 - 1000; 0068 sm.ind1 = find(sm.cond1<=0); 0069 sm.ind2 = find(sm.cond2<=0); 0070 sm.cond1 = 100*(sm.cond1./sm.hoffsetvol); 0071 sm.cond2 = 100*(sm.cond2./sm.hoffsetvol); 0072 sm.sm = (sm.cond1 - sm.cond2)./(sm.cond1 + sm.cond2); 0073 case 4, 0074 sm.cond1 = sm.cond1 - 1000; 0075 sm.cond2 = sm.cond2 - 1000; 0076 sm.cond1 = 100*(sm.cond1./sm.hoffsetvol); 0077 sm.cond2 = 100*(sm.cond2./sm.hoffsetvol); 0078 sm.sm = abs(sm.cond1 - sm.cond2)./(sm.cond1); 0079 case 5, 0080 sm.cond1 = sm.cond1 - 1000; 0081 sm.cond2 = sm.cond2 - 1000; 0082 0083 indless = find(sm.cond1 < sm.cond2); 0084 0085 maxmat = max(sm.cond1, sm.cond2); 0086 minmat = min(sm.cond1, sm.cond2); 0087 0088 sm.sm = (maxmat - minmat)./(maxmat); 0089 sm.sm(indless) = (-1)*sm.sm(indless); 0090 case 6, 0091 sm.cond1 = sm.cond1 - 1000; 0092 sm.cond2 = sm.cond2 - 1000; 0093 sm.cond1 = 100*(sm.cond1./sm.hoffsetvol); 0094 sm.cond2 = 100*(sm.cond2./sm.hoffsetvol); 0095 sm.sm = (exp(sm.cond1) - exp(sm.cond2))./(exp(sm.cond1) + exp(sm.cond2)); 0096 end 0097 0098 clear sm.hoffsetvol; 0099 0100 %% either mask by an ROI or by a fraction 0101 if (sm.ROIs.roi) 0102 fprintf('view_sm: reading in the roi named %s...\n', sm.ROIs.roi); 0103 roi_path = sprintf('%s/%s.roi', sm.ROIs.maskDir, sm.ROIs.roi); 0104 roi = textread(roi_path); 0105 roi = expandROI(roi, size(sm.base)); 0106 0107 sm.virginsm = sm.sm; 0108 sm.sm = sm.sm .* roi; 0109 0110 ind = find(sm.sm ~= 0); 0111 indrest = setdiff(find(sm.base), ind); 0112 0113 sm.sm(indrest) = 0; 0114 sm.base(ind) = 0; 0115 0116 if (sm.ROIs.exclude) 0117 %% find mean and standard deviation of selectivity map 0118 %% we will use this to exclude indices that are x number 0119 %% of standard deviations away from the mean 0120 sm.mean = mean(reshape1d(sm.sm(ind))); 0121 sm.std = std(reshape1d(sm.sm(ind))); 0122 0123 indexclude = find((sm.sm <= (sm.mean - 2*sm.std)) | (sm.sm >= (sm.mean + 2*sm.std))); 0124 0125 sm.sm(indexclude) = 0; 0126 end 0127 0128 %% find the min and the max now, after we have (perhaps) 0129 %% excluded the outlier voxels 0130 sm.min = min(min(min(sm.sm))); 0131 sm.max = max(max(max(sm.sm))); 0132 else 0133 %% find the indices for the values that are both inside and 0134 %% outside our range 0135 sm.frac = 1 - sm.ROIs.frac; 0136 sm.virginsm = sm.sm; 0137 0138 sm.min = min(min(min(sm.sm))); 0139 sm.max = max(max(max(sm.sm))); 0140 if (sm.ROIs.min & sm.ROIs.max) 0141 indmin = find(sm.sm <= sm.ROIs.min); 0142 indmax = find(sm.sm >= sm.ROIs.max); 0143 else 0144 indmin = find(sm.sm <= (sm.min * sm.frac)); 0145 indmax = find(sm.sm >= (sm.max * sm.frac)); 0146 end 0147 0148 ind = [indmin; indmax]; 0149 indrest = setdiff(find(sm.sm), ind); 0150 0151 %% set the values of the selectivity map outside the range to 0, 0152 %% while set the values of the base inside the range to 0. 0153 %% this is necessary for the colormap trickery employed below. 0154 sm.sm(indrest) = 0; 0155 sm.base(ind) = 0; 0156 end 0157 0158 %% store these matrices as originals 0159 sm.originalsm = sm.sm; 0160 sm.originalvirginsm = sm.virginsm; 0161 sm.originalbase = sm.base; 0162 0163 %% error checking 0164 if (~strncmp(sm.display,'mosaic',6) & ~(sm.display <= sm.ROIs.numSlices)) 0165 oops(254); 0166 end 0167 0168 %% create a colormap 128 elements long 0169 %% first 64, the nmr colormap with positive and negative tails 0170 %% second 64, a grayscale colormap which will be used for the base 0171 sm.ncmap = 64; 0172 %%sm.cmap = nmrcolormap(sm.ncmap, 'both'); 0173 %%sm.cmap = jet(sm.ncmap); 0174 sm.cmap = autumn(sm.ncmap); 0175 sm.cmap = [sm.cmap; gray(sm.ncmap)]; 0176 0177 switch sm.display 0178 case 'mosaic', 0179 %% start from scratch 0180 sm.virginsm = sm.originalvirginsm; 0181 sm.sm = sm.originalsm; 0182 sm.base = sm.originalbase; 0183 0184 %% scale the base 0185 %% this scales it so that it fits into the upper part of 0186 %% the colormap 0187 sm.range = range([sm.min sm.max]); 0188 sm.base = rescaleVolume(sm.base, (maxminfactor*sm.max), ((maxminfactor*sm.max)+sm.range)); 0189 0190 %% set the base voxels that will contain a sm value to 0 0191 sm.base(ind) = 0; 0192 0193 %% reshape the data into mosaics 0194 sm.virginsm = vol2mos(permute(sm.virginsm, [2 3 1])); 0195 sm.base = vol2mos(permute(sm.base, [2 3 1])); 0196 sm.sm = vol2mos(permute(sm.sm, [2 3 1])); 0197 0198 %% combine base and overlapROI together 0199 sm.composite = sm.base + sm.sm; 0200 0201 %% display result 0202 sm.himage = image(sm.composite, 'CDataMapping', 'scaled'); 0203 caxis([sm.min*0.99 sm.max+sm.range]); %% set the axis to 0204 %% use the range employed 0205 colorbar('vert'); 0206 sm.hfig = gcf; 0207 sm.haxis = gca; 0208 colormap(sm.cmap); 0209 set(gcf,'UserData',sm); 0210 return; 0211 otherwise, 0212 %% start from scratch 0213 sm.virginsm = sm.originalvirginsm; 0214 sm.sm = sm.originalsm; 0215 sm.base = sm.originalbase; 0216 0217 %% scale the base 0218 %% this scales it so that it fits into the upper part of 0219 %% the colormap 0220 sm.range = range([sm.min sm.max]); 0221 sm.base = rescaleVolume(sm.base, (maxminfactor*sm.max), ((maxminfactor*sm.max)+sm.range)); 0222 0223 %% set the base voxels that will contain a sm value to 0 0224 sm.base(ind) = 0; 0225 0226 %% get the correct slice 0227 sm.virginsm = squeeze(sm.virginsm(sm.display,:,:)); 0228 sm.base = squeeze(sm.base(sm.display,:,:)); 0229 sm.sm = squeeze(sm.sm(sm.display,:,:)); 0230 0231 %% combine base and overlapROI together 0232 sm.composite = sm.base + sm.sm; 0233 0234 %% display result 0235 sm.himage = image(sm.composite, 'CDataMapping', 'scaled'); 0236 caxis([sm.min*0.99 sm.max+sm.range]); %% set the axis to 0237 %% use the range employed 0238 colorbar('vert'); 0239 sm.hfig = gcf; 0240 sm.haxis = gca; 0241 colormap(sm.cmap); 0242 set(gcf,'UserData',sm); 0243 return; 0244 end 0245 end 0246 0247 %% parse the callback function 0248 sm = get(gcf,'UserData'); 0249 flag = deblank(varargin{1}); 0250 0251 switch(flag) 0252 case {'kbd'}, 0253 c = get(sm.hfig,'CurrentCharacter'); 0254 switch(c) 0255 case 's', %% want to change a slice 0256 title = 'Enter a slice number'; 0257 prompt = {'Slice number:'}; 0258 lines = 1; 0259 def = {'1'}; 0260 answer = inputdlg(prompt,title,lines,def); 0261 0262 %% if we actually have a slice number 0263 if (~isempty(answer)) 0264 sliceNum = sscanf(answer{1},'%d'); 0265 if ((sliceNum <= 0) | (sliceNum > sm.ROIs.numSlices)) 0266 msg = sprintf('Please enter a valid slice'); 0267 errordlg(msg); 0268 end 0269 fprintf('setting display to slice %d\n',sliceNum); 0270 fprintf('\n'); 0271 view_sm(sm.ROIs,sm.type,sliceNum); 0272 end 0273 case 'm', %% set to mosaic 0274 fprintf('setting display to mosaic\n'); 0275 fprintf('\n'); 0276 view_sm(sm.ROIs,sm.type,'mosaic'); 0277 case 'h', %% give help info 0278 s = sprintf('Keypress Commands:\n'); 0279 s = sprintf('%sh - this help information\n',s); 0280 s = sprintf('%sm - display mosaic view\n',s); 0281 s = sprintf('%sq - quit\n',s); 0282 s = sprintf('%ss - select a specific slice\n',s); 0283 s = sprintf('%s\nTo see a selectivity index\n',s); 0284 s = sprintf('%sfor a particular voxel, click on\n',s); 0285 s = sprintf('%sthe voxel. If the selecivity\n',s); 0286 s = sprintf('%svalue for the voxel lies outside\n',s); 0287 s = sprintf('%sthe range specified on the command\n',s); 0288 s = sprintf('%sline, no selectivity index will\n',s); 0289 s = sprintf('%sbe given.',s); 0290 msgbox(s,'view_sm Help','help'); 0291 fprintf('\n'); 0292 case 'q', %% quit 0293 hoverlap = gcf; 0294 close(hoverlap); 0295 return; 0296 end 0297 0298 case {'wbd'}, 0299 xyz = get(gca, 'CurrentPoint'); 0300 smval = sm.virginsm(round(xyz(1,2)), round(xyz(1,1))); 0301 %%if (smval == 0) 0302 %% fprintf('Selectivity index does not fall within the range specified on the command line.'); 0303 %%else 0304 fprintf('Selectivity Index: %f', smval); 0305 %%end 0306 fprintf('\n'); 0307 end