[80] | 1 | function map_struct = aedes_fitmaps(DATA,maptype,fit_vals,varargin) |
---|
| 2 | % AEDES_FITMAPS - Calculate various parameter maps (T1, T2, T1rho, B1, |
---|
[57] | 3 | % Perfusion, etc.) |
---|
[2] | 4 | % |
---|
| 5 | % |
---|
| 6 | % Synopsis: |
---|
[80] | 7 | % MapStruct = aedes_fitmaps(data,maptype,fit_vals,... |
---|
[29] | 8 | % 'property1',value1,'property2',value2,...) |
---|
[2] | 9 | % |
---|
| 10 | % Description: |
---|
[57] | 11 | % Function fits various parameter maps defined by MAPTYPE from the |
---|
| 12 | % slice data DATA using the fit values FIT_VALS. The DATA can be a 3D |
---|
| 13 | % matrix or an Aedes DATA-structure. The MAPTYPE variable is a string |
---|
| 14 | % variable that defines the type of the fitted map, i.e. the function |
---|
| 15 | % used in the fitting. Valid values for MAPTYPE are the following: |
---|
[2] | 16 | % |
---|
[29] | 17 | % 'T1_IR' <-> T1 map (inversion recovery) |
---|
| 18 | % 'T1_SR' <-> T1 map (saturation recovery) |
---|
| 19 | % 'T1_3P' <-> T1 map (3-parameter fit) |
---|
[2] | 20 | % 'T2' <-> T2 map |
---|
[18] | 21 | % 'R2' <-> R2 map |
---|
[2] | 22 | % 'T1r' <-> T1 rho map |
---|
[30] | 23 | % 'T2r' <-> T2 rho map |
---|
| 24 | % 'ADC' <-> Apparent diffusion coefficient map |
---|
[2] | 25 | % 'perfusion' <-> Perfusion map |
---|
| 26 | % 'B1' <-> B1 map |
---|
| 27 | % |
---|
[24] | 28 | % If you want to omit some slices from the map calculations, just |
---|
| 29 | % set the corresponding value in FIT_VALS to NaN. |
---|
| 30 | % |
---|
[29] | 31 | % Property-value pairs can be used to control additional options in |
---|
| 32 | % the fitting of maps (The { } denotes default value for the |
---|
| 33 | % corresponding property): |
---|
| 34 | % |
---|
| 35 | % Property: Value: Description: |
---|
| 36 | % ******** ******** ************ |
---|
| 37 | % 'FileName' String % Save maps to files |
---|
| 38 | % |
---|
[34] | 39 | % 'SaveSpinDensities' ['on' | {'off'}] % Save also the |
---|
| 40 | % % S0-parameter in a |
---|
| 41 | % % separate file. |
---|
| 42 | % % This property is ignored |
---|
| 43 | % % if the FileName -property |
---|
| 44 | % % is not defined |
---|
| 45 | % |
---|
[29] | 46 | % 'Wbar' [{'on'} | 'off' ] % Show/hide waitbar |
---|
| 47 | % |
---|
| 48 | % 'Linear' [{'on'} | 'off'] % Perform linear or |
---|
| 49 | % % non-linear fit |
---|
| 50 | % |
---|
| 51 | % 'InitVal' vector % Initial values for |
---|
| 52 | % % non-linear fit |
---|
| 53 | % |
---|
| 54 | % 'MaxIter' scalar % Maximum number of |
---|
| 55 | % % iterations in non-linear |
---|
[34] | 56 | % % fits (default=50) |
---|
[29] | 57 | % |
---|
| 58 | % The 'FileName' property can be used to write the map into a file |
---|
[51] | 59 | % with the corresponding file extension (.t1, .t2, .t1r, etc.). The |
---|
| 60 | % files are normal Matlab MAT-files with a different file extension |
---|
| 61 | % and all the parameters used to calculate the map are also stored in |
---|
| 62 | % the file. Aedes can read these files normally. You can also load |
---|
| 63 | % these files into Matlab workspace by typing |
---|
| 64 | % maps=load('/path/to/my/mapfile','-mat'). If the 'FileName' property |
---|
[34] | 65 | % does not include full path, the map-files a written into the same |
---|
| 66 | % directory as the data by default if possible. Otherwise the maps |
---|
| 67 | % are written into the current directory. |
---|
[29] | 68 | % |
---|
| 69 | % The 'Linear' property is ignored in the cases where only linear or |
---|
| 70 | % non-linear fit is available. |
---|
| 71 | % |
---|
| 72 | % The 'InitVal' property contains the global initial values that are |
---|
| 73 | % used for all individual fittings. If the 'InitVal' property is |
---|
| 74 | % omitted, some "optimal" initial values are estimated using the data |
---|
| 75 | % and the fit values. |
---|
| 76 | % |
---|
[2] | 77 | % The function outputs a structure containing the following fields: |
---|
| 78 | % |
---|
| 79 | % MapStruct |
---|
| 80 | % |-> Map :(The calculated map(s)) |
---|
| 81 | % |-> S0 :(The "spin density") |
---|
| 82 | % |-> FitValues :(Values used in the fitting) |
---|
| 83 | % |-> Type :(Used map type, e.g. 'T2') |
---|
| 84 | % |-> Linear :(1=linear fitting, 0=nonlinear fitting) |
---|
[30] | 85 | % |-> MaxIter :(Maximum number of iterations, non-linear only) |
---|
| 86 | % |-> ParentFileName :(File(s) from which the map was calculated) |
---|
[2] | 87 | % |
---|
| 88 | % |
---|
| 89 | % Examples: |
---|
| 90 | % |
---|
| 91 | % See also: |
---|
[37] | 92 | % AEDES |
---|
[2] | 93 | |
---|
[37] | 94 | % This function is a part of Aedes - A graphical tool for analyzing |
---|
[36] | 95 | % medical images |
---|
[2] | 96 | % |
---|
[36] | 97 | % Copyright (C) 2006 Juha-Pekka Niskanen <Juha-Pekka.Niskanen@uku.fi> |
---|
| 98 | % |
---|
[45] | 99 | % Department of Physics, Department of Neurobiology |
---|
[39] | 100 | % University of Kuopio, FINLAND |
---|
[36] | 101 | % |
---|
| 102 | % This program may be used under the terms of the GNU General Public |
---|
| 103 | % License version 2.0 as published by the Free Software Foundation |
---|
| 104 | % and appearing in the file LICENSE.TXT included in the packaging of |
---|
| 105 | % this program. |
---|
| 106 | % |
---|
| 107 | % This program is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE |
---|
| 108 | % WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. |
---|
[2] | 109 | |
---|
| 110 | map_struct=[]; |
---|
| 111 | if nargin<3 |
---|
| 112 | error('Too few input arguments!') |
---|
| 113 | end |
---|
| 114 | |
---|
[30] | 115 | Dat.filename = ''; |
---|
| 116 | Dat.linear = true; |
---|
| 117 | Dat.wbar = true; |
---|
[129] | 118 | Dat.MaxIter = 75; % Maximum number of iterations in nonlinear fits |
---|
[30] | 119 | Dat.initVal = []; % Empty means that defaults are used for initial values |
---|
| 120 | ParentFileName = ''; |
---|
| 121 | Dat.UseComplexData = false; |
---|
[34] | 122 | Dat.SaveSpinDensities = false; |
---|
[73] | 123 | Dat.Mask = []; |
---|
[2] | 124 | |
---|
[29] | 125 | % Parse additional input arguments |
---|
| 126 | for ii=1:2:length(varargin) |
---|
| 127 | switch lower(varargin{ii}) |
---|
[73] | 128 | case 'filename' |
---|
| 129 | Dat.filename = varargin{ii+1}; |
---|
| 130 | case 'linear' |
---|
| 131 | if ischar(varargin{ii+1}) |
---|
| 132 | if strcmpi(varargin{ii+1},'off') |
---|
| 133 | Dat.linear = false; |
---|
| 134 | end |
---|
| 135 | elseif islogical(varargin{ii+1}) |
---|
| 136 | Dat.linear = varargin{ii+1}; |
---|
| 137 | else |
---|
| 138 | error('Invalid value for the LINEAR property!'); |
---|
| 139 | end |
---|
| 140 | case {'wbar','waitbar'} |
---|
| 141 | if ischar(varargin{ii+1}) |
---|
| 142 | if strcmpi(varargin{ii+1},'off') |
---|
| 143 | Dat.wbar = false; |
---|
| 144 | end |
---|
| 145 | elseif islogical(varargin{ii+1}) |
---|
| 146 | Dat.wbar = varargin{ii+1}; |
---|
| 147 | else |
---|
| 148 | error('Invalid value for the WBAR property!'); |
---|
| 149 | end |
---|
| 150 | case {'initval','initialvalues'} |
---|
| 151 | Dat.initVal = varargin{ii+1}; |
---|
| 152 | case 'maxiter' |
---|
| 153 | Dat.MaxIter = varargin{ii+1}; |
---|
| 154 | case 'usecomplexdata' |
---|
| 155 | Dat.UseComplexData = varargin{ii+1}; |
---|
| 156 | case 'savespindensities' |
---|
| 157 | if ischar(varargin{ii+1}) |
---|
| 158 | if strcmpi(varargin{ii+1},'on') |
---|
| 159 | Dat.SaveSpinDensities = true; |
---|
| 160 | end |
---|
| 161 | elseif islogical(varargin{ii+1}) |
---|
| 162 | Dat.SaveSpinDensities = varargin{ii+1}; |
---|
| 163 | else |
---|
| 164 | error('Invalid value for the WBAR property!'); |
---|
| 165 | end |
---|
| 166 | case 'mask' |
---|
| 167 | Dat.Mask = varargin{ii+1}; |
---|
| 168 | otherwise |
---|
| 169 | error('Invalid property "%s"!',varargin{ii}) |
---|
[2] | 170 | end |
---|
| 171 | end |
---|
| 172 | |
---|
| 173 | if iscell(DATA) && length(DATA)==1 |
---|
| 174 | DATA=DATA{1}; |
---|
| 175 | end |
---|
| 176 | |
---|
| 177 | %% Check map type |
---|
| 178 | if ~ischar(maptype) |
---|
| 179 | error(['Map type has to be one of the following strings: ',... |
---|
[31] | 180 | '%s, %s, %s, %s, %s, %s, %s.'],... |
---|
| 181 | 'T1','T2','T1r','T2r','perfusion','ADC','B1') |
---|
[2] | 182 | end |
---|
| 183 | |
---|
| 184 | %% Convert data to 3D-matrix |
---|
| 185 | if iscell(DATA) && length(DATA)>1 |
---|
| 186 | sz = size(DATA{1}.FTDATA); |
---|
| 187 | sz(3) = length(DATA); |
---|
| 188 | datamtx = zeros(sz); |
---|
[30] | 189 | ParentFileName = {}; |
---|
[2] | 190 | for ii=1:sz(3) |
---|
| 191 | datamtx(:,:,ii)=double(DATA{ii}.FTDATA); |
---|
[30] | 192 | ParentFileName{ii} = fullfile(DATA{ii}.HDR.fpath,DATA{ii}.HDR.fname); |
---|
[2] | 193 | end |
---|
| 194 | elseif isstruct(DATA) |
---|
[62] | 195 | if ndims(DATA.FTDATA)==4 |
---|
| 196 | DATA.FTDATA = squeeze(DATA.FTDATA); |
---|
| 197 | end |
---|
[2] | 198 | sz = size(DATA.FTDATA); |
---|
| 199 | datamtx = double(DATA.FTDATA); |
---|
[30] | 200 | ParentFileName = fullfile(DATA.HDR.fpath,DATA.HDR.fname); |
---|
[2] | 201 | else |
---|
[62] | 202 | if ndims(DATA)==4 |
---|
| 203 | DATA = squeeze(DATA); |
---|
| 204 | end |
---|
[2] | 205 | sz = size(DATA); |
---|
| 206 | datamtx = double(DATA); |
---|
| 207 | end |
---|
| 208 | |
---|
| 209 | %% Check the number of maps to be calculated |
---|
| 210 | if rem(sz(3),length(fit_vals))~=0 |
---|
| 211 | error('Number of slices doesn''t match with the length of fit values.'); |
---|
| 212 | return |
---|
[42] | 213 | elseif length(fit_vals)<2 |
---|
| 214 | % Display error if the number of fit values is not acceptable |
---|
| 215 | error('Cannot calculate maps with less than 2 fit values!') |
---|
[2] | 216 | else |
---|
[30] | 217 | Dat.nMaps = sz(3)/length(fit_vals); |
---|
[24] | 218 | |
---|
| 219 | %% NaN:s in fit_vals means that the corresponding slices will be omitted |
---|
| 220 | %% from the calculations |
---|
| 221 | indNans = find(isnan(fit_vals)); |
---|
| 222 | if not(isempty(indNans)) |
---|
[73] | 223 | l=length(fit_vals); |
---|
| 224 | ind = repmat(indNans,1,Dat.nMaps)+reshape(repmat(l.*[0:Dat.nMaps-1],length(indNans),1),[],1)'; |
---|
| 225 | datamtx(:,:,ind) = []; |
---|
| 226 | fit_vals(indNans)=[]; |
---|
| 227 | sz=size(datamtx); |
---|
[24] | 228 | end |
---|
| 229 | |
---|
[2] | 230 | end |
---|
| 231 | |
---|
[73] | 232 | %% Check that mask is of correct size |
---|
| 233 | if ~isempty(Dat.Mask) |
---|
| 234 | if ~( size(Dat.Mask,1)==sz(1) && ... |
---|
| 235 | size(Dat.Mask,2)==sz(2) && ... |
---|
| 236 | size(Dat.Mask,3)==Dat.nMaps) |
---|
| 237 | error('Mask size does not correspond with data size!') |
---|
| 238 | end |
---|
| 239 | else |
---|
| 240 | Dat.Mask = true(sz(1),sz(2),Dat.nMaps); |
---|
| 241 | end |
---|
[2] | 242 | |
---|
| 243 | switch lower(maptype) |
---|
| 244 | %% Calculate ADC-maps ---------------------------- |
---|
[31] | 245 | case {'adc','df','diffusion'} |
---|
[2] | 246 | |
---|
| 247 | % Fit ADC map |
---|
[30] | 248 | map_out=l_ADC_map(datamtx,fit_vals,'',Dat); |
---|
[2] | 249 | map_out.Type = 'ADC'; |
---|
| 250 | map_out.Linear = true; |
---|
[31] | 251 | fext{1} = 'df'; |
---|
| 252 | fext{2} = 'sf'; |
---|
[2] | 253 | |
---|
| 254 | %% Calculate T1-maps ---------------------------- |
---|
[29] | 255 | case {'t1_ir','t1_sr','t1_3p','t1ir','t1sr','t13p'} |
---|
[2] | 256 | |
---|
[29] | 257 | if any(strcmpi(maptype,{'t1_ir','t1ir'})) |
---|
| 258 | maptype = 'T1_IR'; |
---|
| 259 | elseif any(strcmpi(maptype,{'t1_sr','t1sr'})) |
---|
| 260 | maptype = 'T1_SR'; |
---|
| 261 | else |
---|
| 262 | maptype = 'T1_3P'; |
---|
| 263 | end |
---|
[24] | 264 | |
---|
| 265 | % Fit T1 maps |
---|
[30] | 266 | map_out = l_T1_map(datamtx,fit_vals,maptype,Dat); |
---|
[29] | 267 | map_out.Type = maptype; |
---|
[30] | 268 | map_out.Linear = false;%Dat.linear; |
---|
| 269 | map_out.MaxIter = Dat.MaxIter; |
---|
[24] | 270 | fext{1} = 't1'; |
---|
| 271 | fext{2} = 's1'; |
---|
[2] | 272 | |
---|
| 273 | %% Calculate T1rho-maps ---------------------------- |
---|
| 274 | case {'t1r','t1rho'} |
---|
| 275 | |
---|
[30] | 276 | maptype = 'T1r'; |
---|
| 277 | |
---|
| 278 | % Fit T1r maps |
---|
[34] | 279 | map_out = l_T2_map(datamtx,fit_vals,maptype,Dat); |
---|
| 280 | map_out.Type = 'T1r'; |
---|
[30] | 281 | map_out.Linear = Dat.linear; |
---|
[2] | 282 | fext{1} = 't1r'; |
---|
| 283 | fext{2} = 's1r'; |
---|
| 284 | |
---|
[30] | 285 | %% Calculate T2rho-maps ---------------------------- |
---|
| 286 | case {'t2r','t2rho'} |
---|
[2] | 287 | |
---|
[30] | 288 | maptype = 'T2r'; |
---|
| 289 | |
---|
[34] | 290 | % Fit T2r maps |
---|
| 291 | map_out = l_T2_map(datamtx,fit_vals,maptype,Dat); |
---|
| 292 | map_out.Type = 'T2r'; |
---|
| 293 | map_out.Linear = Dat.linear; |
---|
| 294 | fext{1} = 't2r'; |
---|
| 295 | fext{2} = 's2r'; |
---|
[30] | 296 | |
---|
[34] | 297 | |
---|
[2] | 298 | %% Calculate T2 maps ---------------------------- |
---|
[18] | 299 | case {'t2','r2'} |
---|
[2] | 300 | |
---|
[18] | 301 | if strcmpi(maptype,'t2') |
---|
| 302 | % Fit T2 maps |
---|
[34] | 303 | map_out = l_T2_map(datamtx,fit_vals,maptype,Dat); |
---|
[18] | 304 | map_out.Type = 'T2'; |
---|
| 305 | map_out.Linear = true; |
---|
| 306 | fext{1} = 't2'; |
---|
| 307 | fext{2} = 's2'; |
---|
| 308 | elseif strcmpi(maptype,'r2') |
---|
| 309 | % Fit R2 maps |
---|
[34] | 310 | map_out = l_T2__map(datamtx,fit_vals,maptype,Dat); |
---|
[18] | 311 | map_out.Type = 'R2'; |
---|
[30] | 312 | map_out.Linear = Dat.linear; |
---|
[18] | 313 | map_out.Map = 1./map_out.Map; |
---|
| 314 | map_out.Map(isinf(map_out.Map))=0; |
---|
| 315 | map_out.Map(isnan(map_out.Map))=0; |
---|
| 316 | fext{1} = 'r2'; |
---|
| 317 | fext{2} = 's2'; |
---|
| 318 | end |
---|
[2] | 319 | |
---|
| 320 | |
---|
| 321 | %% Calculate perfusion-maps ---------------------------- |
---|
| 322 | case {'perfusion','perf'} |
---|
| 323 | |
---|
| 324 | |
---|
| 325 | |
---|
| 326 | %% Calculate B1-maps ---------------------------- |
---|
| 327 | case 'b1' |
---|
[30] | 328 | maptype = 'B1'; |
---|
| 329 | |
---|
| 330 | % Fit T1r maps |
---|
[55] | 331 | if isstruct(DATA) && isfield(DATA,'KSPACE') && ... |
---|
| 332 | ~isempty(DATA.KSPACE) && Dat.UseComplexData |
---|
| 333 | datamtx = double(DATA.KSPACE); |
---|
| 334 | datamtx=fftshift(fftshift(fft(fft(datamtx,[],1),[],2),1),2); |
---|
| 335 | fprintf(1,'Calculating B1-maps using complex data...\n') |
---|
| 336 | map_out = l_B1_map(datamtx,fit_vals,maptype,Dat); |
---|
| 337 | else |
---|
| 338 | map_out = l_B1_map(datamtx,fit_vals,maptype,Dat); |
---|
| 339 | end |
---|
[30] | 340 | map_out.Type = maptype; |
---|
| 341 | map_out.Linear = false;%Dat.linear; |
---|
| 342 | map_out.MaxIter = Dat.MaxIter; |
---|
[55] | 343 | fext{1} = 'b1'; |
---|
[30] | 344 | fext{2} = 'mat'; |
---|
[2] | 345 | |
---|
| 346 | case 'mt' |
---|
| 347 | |
---|
[30] | 348 | map_out=l_MTmap(datamtx,fit_vals,Dat); |
---|
[2] | 349 | S0=[]; |
---|
| 350 | |
---|
| 351 | otherwise |
---|
| 352 | error('Unknown map type "%s"!',maptype) |
---|
| 353 | end |
---|
| 354 | |
---|
[30] | 355 | % Add name of the parent file to the structure |
---|
| 356 | map_out.ParentFileName = ParentFileName; |
---|
[29] | 357 | |
---|
[2] | 358 | %% Write maps to files |
---|
[30] | 359 | if ~isempty(Dat.filename) |
---|
[29] | 360 | if nargout==0 |
---|
| 361 | clear map_struct |
---|
| 362 | end |
---|
| 363 | |
---|
[30] | 364 | [fp,fn,fe]=fileparts(Dat.filename); |
---|
[29] | 365 | if isempty(fp) |
---|
[34] | 366 | if ~isempty(ParentFileName) |
---|
| 367 | if iscell(ParentFileName) |
---|
| 368 | [p,n,e]=fileparts(ParentFileName{1}); |
---|
| 369 | fpath = [p,filesep]; |
---|
| 370 | else |
---|
| 371 | [p,n,e]=fileparts(ParentFileName); |
---|
| 372 | fpath = [p,filesep]; |
---|
| 373 | end |
---|
| 374 | else |
---|
| 375 | fpath = [pwd,filesep]; |
---|
| 376 | end |
---|
[29] | 377 | else |
---|
| 378 | fpath = [fp,filesep]; |
---|
| 379 | end |
---|
| 380 | if isempty(fn) |
---|
| 381 | fname = 'mapdata'; |
---|
| 382 | else |
---|
| 383 | fname = fn; |
---|
| 384 | end |
---|
[2] | 385 | for ii=1:size(map_out.Map,3) |
---|
| 386 | Data = map_out.Map(:,:,ii); |
---|
[47] | 387 | Param = map_out; |
---|
| 388 | Param = rmfield(Param,'Map'); |
---|
[2] | 389 | |
---|
| 390 | % Save map |
---|
| 391 | filename = sprintf('%s%s_%03d.%s',fpath,fname,ii,fext{1}); |
---|
| 392 | save(filename,'Data','Param','-mat') |
---|
| 393 | |
---|
| 394 | % Save "spin densities" |
---|
[34] | 395 | if isfield(map_out,'S0') && Dat.SaveSpinDensities |
---|
[30] | 396 | Data = map_out.S0(:,:,ii); |
---|
| 397 | filename = sprintf('%s%s_%03d.%s',fpath,fname,ii,fext{2}); |
---|
| 398 | save(filename,'Data','Param','-mat') |
---|
| 399 | end |
---|
[2] | 400 | end |
---|
| 401 | else |
---|
| 402 | map_struct = map_out; |
---|
| 403 | end |
---|
| 404 | |
---|
| 405 | |
---|
| 406 | |
---|
| 407 | %%%%%%%%%%%%%%%%%%%%% |
---|
| 408 | % Calculate T1-map |
---|
| 409 | %%%%%%%%%%%%%%%%%%%%% |
---|
[30] | 410 | function map_out=l_T1_map(data,TI,maptype,Dat) |
---|
[24] | 411 | % Calculate T1 relaxation times |
---|
| 412 | % |
---|
| 413 | % T1 functions: |
---|
[29] | 414 | % Inversion recovery : S=abs(S0*(1-2*exp(-TI/T1))) |
---|
[24] | 415 | % Saturation recovery : S=S0*(1-1*exp(-TI/T1)) |
---|
| 416 | % 3 Parameter fit : S=S0*(1-A*exp(-TI/T1)) |
---|
| 417 | % |
---|
| 418 | % where: |
---|
| 419 | % S = measured signal, S0 = "spin density" |
---|
| 420 | % TI = inversion time, T1 = T1 relaxation time |
---|
[2] | 421 | |
---|
[29] | 422 | %% Data size |
---|
| 423 | sz=size(data); |
---|
| 424 | |
---|
[30] | 425 | %% Fit values |
---|
[29] | 426 | TI = TI(:); |
---|
| 427 | |
---|
| 428 | %% Data slice index matrix |
---|
[30] | 429 | IndMtx = reshape(1:sz(3),[length(TI) Dat.nMaps])'; |
---|
[29] | 430 | |
---|
[30] | 431 | %% Allocate space for parameters |
---|
| 432 | map = zeros([sz(1) sz(2) Dat.nMaps]); |
---|
| 433 | S0 = zeros([sz(1) sz(2) Dat.nMaps]); |
---|
| 434 | A = zeros([sz(1) sz(2) Dat.nMaps]); |
---|
[29] | 435 | |
---|
[30] | 436 | % Fit options for fminsearch |
---|
[29] | 437 | options = optimset('Display','off',... |
---|
[30] | 438 | 'MaxIter',Dat.MaxIter); |
---|
[29] | 439 | |
---|
| 440 | estimateInitVal = false; |
---|
[30] | 441 | if isempty(Dat.initVal) |
---|
[29] | 442 | estimateInitVal = true; |
---|
[24] | 443 | end |
---|
[2] | 444 | |
---|
[29] | 445 | % - Inversion recovery |
---|
[30] | 446 | if strcmpi(maptype,'t1_ir') |
---|
[29] | 447 | t1_map_type = 1; |
---|
[30] | 448 | elseif strcmpi(maptype,'t1_sr') |
---|
[29] | 449 | % Saturation recovery |
---|
| 450 | t1_map_type = 2; |
---|
| 451 | else |
---|
| 452 | % 3 Parameter fit |
---|
| 453 | t1_map_type = 3; |
---|
| 454 | end |
---|
| 455 | |
---|
[30] | 456 | % Variables for waitbar |
---|
| 457 | nFits = (sz(1)*sz(2)*Dat.nMaps); |
---|
| 458 | counter = 1; |
---|
[51] | 459 | meanTI = mean(TI); |
---|
[29] | 460 | |
---|
[30] | 461 | % Calculate maps |
---|
| 462 | for ii=1:Dat.nMaps |
---|
| 463 | if Dat.wbar && ii==1 |
---|
[80] | 464 | wbh=aedes_wbar(0,sprintf('Processing map %d/%d',ii,Dat.nMaps)); |
---|
[30] | 465 | elseif Dat.wbar |
---|
[80] | 466 | aedes_wbar(counter/nFits,wbh,sprintf('Processing map %d/%d',ii,Dat.nMaps)); |
---|
[29] | 467 | end |
---|
[73] | 468 | |
---|
[29] | 469 | for kk=1:sz(1) |
---|
[73] | 470 | for tt=1:sz(2) |
---|
| 471 | if ~isempty(Dat.Mask) && ~Dat.Mask(kk,tt,ii) |
---|
| 472 | if Dat.wbar |
---|
[80] | 473 | aedes_wbar(counter/nFits,wbh); |
---|
[73] | 474 | end |
---|
| 475 | counter = counter+1; |
---|
| 476 | continue |
---|
| 477 | end |
---|
| 478 | |
---|
| 479 | % Data values |
---|
| 480 | data_val = squeeze(data(kk,tt,IndMtx(ii,:))); |
---|
| 481 | data_val = data_val(:); |
---|
| 482 | |
---|
| 483 | % Initial values for the fit |
---|
| 484 | if estimateInitVal |
---|
| 485 | if t1_map_type == 3 |
---|
| 486 | init_val = [max(data_val); meanTI; 2]; |
---|
| 487 | else |
---|
| 488 | init_val = [max(data_val); meanTI]; |
---|
| 489 | end |
---|
| 490 | else |
---|
| 491 | init_val = Dat.initVal; |
---|
| 492 | end |
---|
| 493 | |
---|
| 494 | % Nelder-Mead simplex iteration |
---|
| 495 | if t1_map_type == 1 |
---|
| 496 | fhandle = @(x) norm(data_val - abs(x(1)*(1-2*exp(-TI./x(2))))); |
---|
| 497 | elseif t1_map_type == 2 |
---|
| 498 | fhandle = @(x) norm(data_val - abs(x(1)*(1-1*exp(-TI./x(2))))); |
---|
| 499 | else |
---|
| 500 | fhandle = @(x) norm(data_val - abs(x(1)*(1-x(3)*exp(-TI./x(2))))); |
---|
| 501 | end |
---|
| 502 | th = fminsearch(fhandle,init_val,options); |
---|
| 503 | |
---|
| 504 | S0(kk,tt,ii) = th(1); |
---|
| 505 | map(kk,tt,ii) = th(2); |
---|
| 506 | if t1_map_type == 3 |
---|
| 507 | A(kk,tt,ii) = th(3); |
---|
| 508 | end |
---|
| 509 | |
---|
| 510 | if Dat.wbar |
---|
[80] | 511 | aedes_wbar(counter/nFits,wbh); |
---|
[73] | 512 | end |
---|
| 513 | counter = counter+1; |
---|
| 514 | end |
---|
[29] | 515 | end |
---|
[24] | 516 | end |
---|
[30] | 517 | if Dat.wbar |
---|
[29] | 518 | close(wbh) |
---|
| 519 | end |
---|
[73] | 520 | |
---|
[29] | 521 | map(map<0)=0; |
---|
| 522 | map(isinf(map))=0; |
---|
| 523 | map(isnan(map))=0; |
---|
[2] | 524 | |
---|
[29] | 525 | map_out.Map = map; |
---|
| 526 | map_out.S0 = S0; |
---|
| 527 | map_out.Angle = A; |
---|
| 528 | map_out.FitValues = TI; |
---|
| 529 | |
---|
| 530 | |
---|
[30] | 531 | |
---|
| 532 | %%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 533 | % Calculate B1-map |
---|
| 534 | %%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 535 | function map_out = l_B1_map(data,fit_vals,maptype,Dat) |
---|
| 536 | |
---|
| 537 | %% Data size |
---|
| 538 | sz=size(data); |
---|
| 539 | |
---|
| 540 | %% Fit values |
---|
| 541 | fit_vals = fit_vals(:); |
---|
| 542 | |
---|
| 543 | %% Data slice index matrix |
---|
| 544 | IndMtx = reshape(1:sz(3),[length(fit_vals) Dat.nMaps])'; |
---|
| 545 | |
---|
| 546 | %% Allocate space for parameters |
---|
| 547 | map = zeros([sz(1) sz(2) Dat.nMaps]); |
---|
| 548 | A = zeros([sz(1) sz(2) Dat.nMaps]); |
---|
| 549 | phase = zeros([sz(1) sz(2) Dat.nMaps]); |
---|
| 550 | |
---|
| 551 | % Fit options for fminsearch |
---|
| 552 | options = optimset('Display','off',... |
---|
| 553 | 'MaxIter',Dat.MaxIter); |
---|
| 554 | |
---|
| 555 | estimateInitVal = false; |
---|
| 556 | if isempty(Dat.initVal) |
---|
| 557 | estimateInitVal = true; |
---|
| 558 | end |
---|
| 559 | |
---|
| 560 | |
---|
| 561 | % Variables for waitbar |
---|
| 562 | nFits = (sz(1)*sz(2)*Dat.nMaps); |
---|
| 563 | counter = 1; |
---|
| 564 | |
---|
| 565 | |
---|
| 566 | % Loop over maps |
---|
| 567 | for ii=1:Dat.nMaps |
---|
| 568 | if Dat.wbar && ii==1 |
---|
[80] | 569 | wbh=aedes_wbar(0,sprintf('Calculating B1-map %d/%d',ii,Dat.nMaps)); |
---|
[30] | 570 | elseif Dat.wbar |
---|
[80] | 571 | aedes_wbar(counter/nFits,wbh,sprintf('Calculating map %d/%d',ii,Dat.nMaps)); |
---|
[30] | 572 | end |
---|
| 573 | |
---|
| 574 | for tt=1:sz(1) |
---|
| 575 | for kk=1:sz(2) |
---|
| 576 | % Data values and initial values for the iteration |
---|
| 577 | if Dat.UseComplexData |
---|
| 578 | data_val = squeeze(data(kk,tt,IndMtx(ii,:))); |
---|
| 579 | data_val = real(data_val); |
---|
| 580 | data_val = data_val(:); |
---|
| 581 | |
---|
| 582 | if estimateInitVal |
---|
| 583 | val1=length(find(diff(data_val<0)<0)); |
---|
| 584 | val2=length(find(diff(data_val<0)>0)); |
---|
| 585 | tot_t = fit_vals(end)-fit_vals(1); |
---|
| 586 | omega = (((val1+val2)/2)/tot_t)*2*pi; |
---|
| 587 | init_val = [2*max(data_val) omega 0.1]; |
---|
| 588 | init_val=init_val(:); |
---|
| 589 | else |
---|
| 590 | init_val = Dat.initVal; |
---|
| 591 | end |
---|
| 592 | |
---|
| 593 | else |
---|
| 594 | data_val = squeeze(data(kk,tt,IndMtx(ii,:))); |
---|
| 595 | data_val = data_val(:); |
---|
| 596 | |
---|
| 597 | if estimateInitVal |
---|
| 598 | val1=length(find(diff(data_val<(max(data_val)*0.5))<0)); |
---|
| 599 | val2=length(find(diff(data_val<(max(data_val)*0.5))>0)); |
---|
| 600 | tot_t = fit_vals(end)-fit_vals(1); |
---|
| 601 | omega = (((val1+val2)/4)/tot_t)*2*pi*0.9; |
---|
| 602 | init_val = [max(data_val)-min(data_val) omega 0.1]; |
---|
| 603 | init_val=init_val(:); |
---|
| 604 | else |
---|
| 605 | init_val = Dat.initVal; |
---|
| 606 | end |
---|
| 607 | end |
---|
| 608 | |
---|
| 609 | % Function handle for fminsearch |
---|
| 610 | if Dat.UseComplexData |
---|
[51] | 611 | fhandle= @(x) norm(data_val - x(1)*cos(x(2)*fit_vals+x(3))); |
---|
[30] | 612 | else |
---|
[51] | 613 | fhandle= @(x) norm(data_val - abs(x(1)*cos(x(2)*fit_vals+x(3)))); |
---|
[30] | 614 | end |
---|
| 615 | th = fminsearch(fhandle,init_val,options); |
---|
| 616 | |
---|
| 617 | map(kk,tt,ii) = abs(th(2))/(2*pi); |
---|
| 618 | A(kk,tt,ii) = th(1); |
---|
| 619 | phase(kk,tt,ii) = th(3); |
---|
| 620 | |
---|
| 621 | if Dat.wbar |
---|
[80] | 622 | aedes_wbar(counter/nFits,wbh); |
---|
[30] | 623 | end |
---|
| 624 | counter = counter+1; |
---|
| 625 | end |
---|
| 626 | end |
---|
| 627 | end |
---|
| 628 | if Dat.wbar |
---|
| 629 | close(wbh) |
---|
| 630 | end |
---|
| 631 | |
---|
| 632 | map(map<0)=0; |
---|
| 633 | map(isinf(map))=0; |
---|
| 634 | map(isnan(map))=0; |
---|
| 635 | |
---|
| 636 | map_out.Map = map; |
---|
| 637 | map_out.A = A; |
---|
| 638 | map_out.phase = phase; |
---|
| 639 | map_out.FitValues = fit_vals; |
---|
| 640 | |
---|
| 641 | |
---|
[34] | 642 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 643 | % Calculate T2-, T1rho- and T2rho-maps |
---|
| 644 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 645 | function map_out=l_T2_map(data,fit_vals,maptype,Dat) |
---|
[2] | 646 | % Calculate T1rho relaxation times |
---|
| 647 | % |
---|
| 648 | % T1rho functions: |
---|
| 649 | % |
---|
| 650 | % Linear: log(S)=log(S0)-SL/T1rho |
---|
| 651 | % Nonlinear: S=S0*exp(-SL/T1rho) |
---|
| 652 | % |
---|
| 653 | % where S=measured signal |
---|
| 654 | % S0 = "spin density" |
---|
| 655 | % SL = spinlock length |
---|
| 656 | % T1rho = T1rho relaxation time |
---|
| 657 | % |
---|
| 658 | % T2-functions: |
---|
| 659 | % |
---|
| 660 | % Linear: log(S)=log(S0)-te/T2 |
---|
| 661 | % Nonlinear: S=S0*exp(-te/T2) |
---|
| 662 | % |
---|
| 663 | % where S=measured signal |
---|
| 664 | % S0 = "spin density" |
---|
| 665 | % te = echo time |
---|
| 666 | % T2 = T2 relaxation time |
---|
| 667 | |
---|
| 668 | done=false; |
---|
| 669 | map=[]; |
---|
| 670 | S0=[]; |
---|
| 671 | |
---|
| 672 | |
---|
| 673 | %% Make sure that fit values are in column vector |
---|
| 674 | fit_vals=fit_vals(:); |
---|
| 675 | |
---|
| 676 | %% Data size |
---|
| 677 | sz=size(data); |
---|
| 678 | |
---|
| 679 | %% Data slice index matrix |
---|
[30] | 680 | IndMtx = reshape(1:sz(3),[length(fit_vals) Dat.nMaps])'; |
---|
[2] | 681 | |
---|
[30] | 682 | map = zeros([sz(1) sz(2) Dat.nMaps]); |
---|
| 683 | S0 = zeros([sz(1) sz(2) Dat.nMaps]); |
---|
| 684 | A = zeros([sz(1) sz(2) Dat.nMaps]); |
---|
[2] | 685 | |
---|
[51] | 686 | % Intensity values should not be less than zero |
---|
| 687 | data(data<0)=0; |
---|
| 688 | |
---|
[2] | 689 | %% Use linearized form (fast) |
---|
[30] | 690 | if Dat.linear |
---|
| 691 | |
---|
| 692 | nFits = sz(2)*Dat.nMaps; |
---|
| 693 | counter = 1; |
---|
| 694 | |
---|
[2] | 695 | H = [ones(size(fit_vals)) -fit_vals]; |
---|
[30] | 696 | for ii=1:Dat.nMaps |
---|
| 697 | if Dat.wbar && ii==1 |
---|
[80] | 698 | wbh=aedes_wbar(0,sprintf('Processing map %d/%d',ii,Dat.nMaps)); |
---|
[30] | 699 | elseif Dat.wbar |
---|
[80] | 700 | aedes_wbar(counter/nFits,wbh,sprintf('Processing map %d/%d',ii,Dat.nMaps)); |
---|
[30] | 701 | end |
---|
| 702 | |
---|
[2] | 703 | for kk=1:sz(2) |
---|
| 704 | tmp=squeeze(data(:,kk,IndMtx(ii,:))).'; |
---|
| 705 | z=log(tmp); |
---|
| 706 | th=H\z; |
---|
| 707 | S0(:,kk,ii)=exp(th(1,:)'); |
---|
| 708 | |
---|
| 709 | % Try to avoid "divide by zero" warnings |
---|
| 710 | |
---|
| 711 | map_tmp = th(2,:)'; |
---|
| 712 | map_tmp(map_tmp==0)=eps; |
---|
| 713 | |
---|
| 714 | map(:,kk,ii)=1./map_tmp; |
---|
[30] | 715 | counter = counter+1; |
---|
| 716 | if Dat.wbar |
---|
[80] | 717 | aedes_wbar(counter/nFits,wbh); |
---|
[30] | 718 | end |
---|
[2] | 719 | end |
---|
| 720 | end |
---|
| 721 | map(map<0)=0; |
---|
| 722 | map(isinf(map))=0; |
---|
| 723 | map(isnan(map))=0; |
---|
| 724 | else |
---|
[30] | 725 | %% Fit nonlinearly using Nelder-Mead Simplex (slow, but more accurate) |
---|
| 726 | estimateInitVal = false; |
---|
| 727 | if isempty(Dat.initVal) |
---|
| 728 | estimateInitVal = true; |
---|
| 729 | end |
---|
| 730 | |
---|
| 731 | % Fit options |
---|
| 732 | options = optimset('Display','off',... |
---|
| 733 | 'MaxIter',Dat.MaxIter); |
---|
| 734 | |
---|
| 735 | nFits = sz(1)*sz(2)*Dat.nMaps; |
---|
| 736 | counter = 1; |
---|
[51] | 737 | meanFV = mean(fit_vals); |
---|
[2] | 738 | |
---|
[30] | 739 | for ii=1:Dat.nMaps |
---|
| 740 | if Dat.wbar && ii==1 |
---|
[80] | 741 | wbh=aedes_wbar(0,sprintf('Processing map %d/%d',ii,Dat.nMaps)); |
---|
[30] | 742 | elseif Dat.wbar |
---|
[80] | 743 | aedes_wbar(counter/nFits,wbh,sprintf('Processing map %d/%d',ii,Dat.nMaps)); |
---|
[30] | 744 | end |
---|
| 745 | |
---|
| 746 | for tt=1:sz(1) |
---|
| 747 | for kk=1:sz(2) |
---|
| 748 | |
---|
| 749 | % Data for the fit |
---|
| 750 | z = squeeze(data(tt,kk,IndMtx(ii,:))); |
---|
| 751 | z=z(:); |
---|
| 752 | |
---|
| 753 | % Initial values for the fit |
---|
| 754 | if estimateInitVal |
---|
| 755 | %init_val = [mean(z); mean(fit_vals); 1]; |
---|
[51] | 756 | init_val = [max(z); meanFV]; |
---|
[30] | 757 | else |
---|
| 758 | init_val = Dat.initVal; |
---|
| 759 | end |
---|
| 760 | |
---|
| 761 | |
---|
| 762 | %fhandle = @(x) sum((z - x(1)*exp(-fit_vals./x(2))+x(3)).^2); |
---|
[51] | 763 | fhandle = @(x) norm(z - x(1)*exp(-fit_vals./x(2))); |
---|
[30] | 764 | th = fminsearch(fhandle,init_val,options); |
---|
| 765 | |
---|
| 766 | S0(tt,kk,ii) = th(1); |
---|
| 767 | map(tt,kk,ii) = th(2); |
---|
| 768 | %A(tt,kk,ii) = th(3); |
---|
| 769 | |
---|
| 770 | if Dat.wbar |
---|
[80] | 771 | aedes_wbar(counter/nFits,wbh); |
---|
[30] | 772 | end |
---|
| 773 | counter = counter+1; |
---|
| 774 | end |
---|
| 775 | end |
---|
| 776 | end |
---|
| 777 | |
---|
[2] | 778 | end |
---|
[30] | 779 | if Dat.wbar |
---|
| 780 | close(wbh) |
---|
| 781 | end |
---|
[2] | 782 | |
---|
[30] | 783 | |
---|
[2] | 784 | map_out.Map = map; |
---|
| 785 | map_out.S0 = S0; |
---|
| 786 | map_out.FitValues = fit_vals; |
---|
| 787 | |
---|
| 788 | %%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 789 | % Calculate Perfusion-map |
---|
| 790 | %%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 791 | function l_Perfusion_map() |
---|
| 792 | |
---|
| 793 | |
---|
| 794 | |
---|
| 795 | %%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 796 | % Calculate Diffusion-map |
---|
| 797 | %%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 798 | function l_Diffusion_map() |
---|
| 799 | |
---|
| 800 | |
---|
| 801 | |
---|
| 802 | %%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 803 | % Calculate ADC-map |
---|
| 804 | %%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
[31] | 805 | function map_out=l_ADC_map(data,b,opt,Dat) |
---|
[2] | 806 | % Calculate Apparent Diffusion Coefficients |
---|
| 807 | % The fitted functions: |
---|
| 808 | % |
---|
| 809 | % Linear: log(S)=log(S0)-b*ADC |
---|
| 810 | % Nonlinear: S=S0*exp(-b*ADC) |
---|
| 811 | % |
---|
| 812 | % where S=measured signal |
---|
| 813 | % S0 = "spin density" |
---|
| 814 | % b = diffusion weighting |
---|
| 815 | % ADC = apparent diffusion coefficient |
---|
| 816 | |
---|
| 817 | if ~exist('opt','var') |
---|
| 818 | % Default to linearized fit |
---|
| 819 | opt = ''; |
---|
| 820 | end |
---|
| 821 | |
---|
| 822 | b=b(:); |
---|
| 823 | |
---|
| 824 | %% Data size |
---|
| 825 | sz=size(data); |
---|
| 826 | |
---|
| 827 | %% Data slice index matrix |
---|
[31] | 828 | IndMtx = reshape(1:sz(3),[length(b) Dat.nMaps])'; |
---|
[2] | 829 | |
---|
[31] | 830 | ADC = zeros([sz(1) sz(2) Dat.nMaps]); |
---|
| 831 | S0 = zeros([sz(1) sz(2) Dat.nMaps]); |
---|
[2] | 832 | |
---|
[51] | 833 | % Intensity values should not be less than zero |
---|
| 834 | data(data<0)=0; |
---|
| 835 | |
---|
[2] | 836 | %% Use linearized form (fast) |
---|
| 837 | if isempty(opt) || strcmpi(opt,'linear') |
---|
| 838 | H = [ones(size(b)) -b]; |
---|
[31] | 839 | for ii=1:Dat.nMaps |
---|
[2] | 840 | for kk=1:sz(2) |
---|
| 841 | tmp=squeeze(data(:,kk,IndMtx(ii,:))).'; |
---|
| 842 | z=log(tmp); |
---|
| 843 | th=H\z; |
---|
| 844 | S0(:,kk,ii)=exp(th(1,:)'); |
---|
| 845 | ADC(:,kk,ii)=th(2,:)'; |
---|
| 846 | end |
---|
| 847 | end |
---|
| 848 | ADC(ADC<0)=0; |
---|
| 849 | ADC(isinf(ADC))=0; |
---|
| 850 | ADC(isnan(ADC))=0; |
---|
| 851 | else |
---|
| 852 | %% Fit using nonlinear LS (slow, but more accurate) |
---|
| 853 | |
---|
| 854 | end |
---|
| 855 | |
---|
| 856 | map_out.Map = ADC; |
---|
| 857 | map_out.S0 = S0; |
---|
[31] | 858 | map_out.FitValues = b; |
---|
[2] | 859 | |
---|
| 860 | |
---|
| 861 | |
---|
| 862 | |
---|
| 863 | |
---|
| 864 | %%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 865 | % Calculate MT-maps |
---|
| 866 | %%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 867 | function [map]=l_MTmap(data,fit_vals,nMaps) |
---|
| 868 | |
---|
| 869 | sz=size(data); |
---|
| 870 | lims = [-700 700]; |
---|
| 871 | |
---|
| 872 | |
---|
| 873 | |
---|
| 874 | [sorted_vals,ind]=sort(fit_vals); |
---|
| 875 | [mx,mx_ind]=max(abs(sorted_vals)); |
---|
| 876 | sorted_vals(mx_ind)=[]; |
---|
| 877 | sorted_vals = sorted_vals(3:end-2); |
---|
| 878 | |
---|
| 879 | |
---|
| 880 | %% Allocate space for map |
---|
| 881 | map=zeros(sz(1),sz(2)); |
---|
| 882 | |
---|
| 883 | lim1 = find(sorted_vals==lims(1)); |
---|
| 884 | lim2 = find(sorted_vals==lims(2)); |
---|
| 885 | |
---|
| 886 | |
---|
| 887 | if false |
---|
| 888 | |
---|
| 889 | for ii=1:sz(1) |
---|
| 890 | for kk=1:sz(2) |
---|
| 891 | vox_data = squeeze(data(ii,kk,:)); |
---|
| 892 | vox_data=vox_data(ind); |
---|
| 893 | vox_data=vox_data./vox_data(mx_ind); |
---|
| 894 | vox_data(mx_ind)=[]; |
---|
| 895 | vox_data = vox_data(3:end-2); |
---|
| 896 | |
---|
| 897 | df=diff([vox_data(lim1) vox_data(lim2)]); |
---|
| 898 | map(ii,kk)=df; |
---|
| 899 | end |
---|
| 900 | end |
---|
| 901 | |
---|
| 902 | else |
---|
| 903 | |
---|
| 904 | zind=find(sorted_vals==0); |
---|
| 905 | fit=sorted_vals(zind-3:zind+3); |
---|
| 906 | fit=fit(:); |
---|
| 907 | %fit(4)=[]; |
---|
| 908 | H = [fit.^2 fit ones(size(fit))]; |
---|
| 909 | |
---|
| 910 | for ii=1:sz(1) |
---|
| 911 | for kk=1:sz(2) |
---|
| 912 | vox_data = squeeze(data(ii,kk,:)); |
---|
| 913 | vox_data=vox_data(ind); |
---|
| 914 | vox_data=vox_data./vox_data(mx_ind); |
---|
| 915 | vox_data(mx_ind)=[]; |
---|
| 916 | vox_data = vox_data(3:end-2); |
---|
| 917 | |
---|
| 918 | %% Do peak picking by fitting a polynomial |
---|
| 919 | z=vox_data(zind-3:zind+3); |
---|
| 920 | %z(4)=[]; |
---|
| 921 | th=H\z; |
---|
| 922 | zz=fit(1):fit(end); |
---|
| 923 | |
---|
| 924 | [mn,mn_ind]=min(polyval(th,zz)); |
---|
| 925 | shift_val=zz(mn_ind); |
---|
| 926 | %map(ii,kk)=shift_val; |
---|
| 927 | %continue |
---|
| 928 | |
---|
| 929 | |
---|
| 930 | %% Do a spline interpolation to the data |
---|
| 931 | %interp_freq = sorted_vals(1):1:sorted_vals(end); |
---|
| 932 | |
---|
| 933 | interp_data = pchip(sorted_vals,vox_data,... |
---|
| 934 | lims-shift_val); |
---|
| 935 | |
---|
| 936 | |
---|
| 937 | %% Shift frequency values |
---|
| 938 | %interp_freq = interp_freq-shift_val; |
---|
| 939 | |
---|
| 940 | % Get indices to limits |
---|
| 941 | %lim1 = find(interp_freq==lims(1)); |
---|
| 942 | %lim2 = find(interp_freq==lims(2)); |
---|
| 943 | |
---|
| 944 | try |
---|
| 945 | df=diff([interp_data]); |
---|
| 946 | catch |
---|
| 947 | df=1; |
---|
| 948 | end |
---|
| 949 | if abs(df)<0.5 |
---|
| 950 | map(ii,kk)=df; |
---|
| 951 | end |
---|
| 952 | |
---|
| 953 | |
---|
| 954 | % $$$ if ii==70 && kk==42 |
---|
| 955 | % $$$ plot(fit,z,'*-',zz,polyval(th,zz),'r') |
---|
| 956 | % $$$ shift_val |
---|
| 957 | % $$$ pause |
---|
| 958 | % $$$ end |
---|
| 959 | % $$$ |
---|
| 960 | % $$$ if abs(interp_freq(mn_ind))<200 |
---|
| 961 | % $$$ |
---|
| 962 | % $$$ % Shift interpolated data |
---|
| 963 | % $$$ interp_freq = interp_freq-interp_freq(mn_ind); |
---|
| 964 | % $$$ |
---|
| 965 | % $$$ % Get indices to limits |
---|
| 966 | % $$$ lim1 = find(interp_freq==lims(1)); |
---|
| 967 | % $$$ lim2 = find(interp_freq==lims(2)); |
---|
| 968 | % $$$ |
---|
| 969 | % $$$ df=diff([interp_data(lim1) interp_data(lim2)]); |
---|
| 970 | % $$$ if abs(df)<0.5 |
---|
| 971 | % $$$ map(ii,kk)=df; |
---|
| 972 | % $$$ end |
---|
| 973 | % $$$ end |
---|
| 974 | end |
---|
| 975 | %disp(num2str(ii)) |
---|
| 976 | end |
---|
| 977 | |
---|
| 978 | end |
---|