source: an2_fitmaps.m @ 63

Last change on this file since 63 was 62, checked in by tjniskan, 11 years ago
  • Added functions for basic fMRI analysis and resting-state
  • Added trend in the "show voxel timeseries" figure. In EPI datas the

the reference image is excluded from the time series and the y-axis
units are changed to "BOLD (%)".

  • Fixed an endiannes bug in an2_readfdf.m
  • Added plugins basic fmri analysis and resting-state analysis

M misclib/fmri_spm_volumes.m
A misclib/fmri_analysis.m
M an2_readfdf.m
M an2_revision.m
M an2_fitmaps.m
M an2_data_read.m
M aedes.m
A plugins/basic_fmri_analysis.m
A plugins/resting_state_fc.m

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

Powered by Trac 1.0.9.Copyright © Juha-Pekka Niskanen 2008