source: aedes_fitmaps.m

Last change on this file was 129, checked in by tjniskan, 7 years ago
  • Added plugin for correlation mapping with filtering/other options
  • Added flipping of 3D ROIs in aedes.m
  • Added possibility of changing ROI color (in the ROI -uimenu)
  • Minor fixes/changes here and there...

M misclib/fmri_spm_volumes.m
M misclib/fmri_filter.m
M aedes_fitmaps.m
M aedes.m
M plugins/save_roi_as_mask.m
M plugins/calc_asl_cbf.m
A plugins/fmri_plugins/correlation_mapping.m
M plugins/fmri_plugins/resting_state_fc.m
M plugins/map_plugins/t1_inversion_recovery.m
M aedes_revision.m

File size: 24.0 KB
Line 
1function map_struct = aedes_fitmaps(DATA,maptype,fit_vals,varargin)
2% AEDES_FITMAPS - Calculate various parameter maps (T1, T2, T1rho, B1,
3%               Perfusion, etc.)
4%   
5%
6% Synopsis:
7%       MapStruct = aedes_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 = 75; % 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;
123Dat.Mask = [];
124
125% Parse additional input arguments
126for ii=1:2:length(varargin)
127  switch lower(varargin{ii})
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})
170  end
171end
172
173if iscell(DATA) && length(DATA)==1
174  DATA=DATA{1};
175end
176
177%% Check map type
178if ~ischar(maptype)
179  error(['Map type has to be one of the following strings: ',...
180        '%s, %s, %s, %s, %s, %s, %s.'],...
181        'T1','T2','T1r','T2r','perfusion','ADC','B1')
182end
183
184%% Convert data to 3D-matrix
185if iscell(DATA) && length(DATA)>1
186  sz = size(DATA{1}.FTDATA);
187  sz(3) = length(DATA);
188  datamtx = zeros(sz);
189  ParentFileName = {};
190  for ii=1:sz(3)
191    datamtx(:,:,ii)=double(DATA{ii}.FTDATA);
192        ParentFileName{ii} = fullfile(DATA{ii}.HDR.fpath,DATA{ii}.HDR.fname);
193  end
194elseif isstruct(DATA)
195  if ndims(DATA.FTDATA)==4
196    DATA.FTDATA = squeeze(DATA.FTDATA);
197  end
198  sz = size(DATA.FTDATA);
199  datamtx = double(DATA.FTDATA);
200  ParentFileName = fullfile(DATA.HDR.fpath,DATA.HDR.fname);
201else
202  if ndims(DATA)==4
203    DATA = squeeze(DATA);
204  end
205  sz = size(DATA);
206  datamtx = double(DATA);
207end
208
209%% Check the number of maps to be calculated
210if rem(sz(3),length(fit_vals))~=0
211  error('Number of slices doesn''t match with the length of fit values.');
212  return
213elseif 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!')
216else
217  Dat.nMaps = sz(3)/length(fit_vals);
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))
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);
228  end
229 
230end
231
232%% Check that mask is of correct size
233if ~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
239else
240  Dat.Mask = true(sz(1),sz(2),Dat.nMaps);
241end
242
243switch lower(maptype)
244  %% Calculate ADC-maps ----------------------------
245  case {'adc','df','diffusion'}
246   
247    % Fit ADC map
248    map_out=l_ADC_map(datamtx,fit_vals,'',Dat);
249    map_out.Type = 'ADC';
250    map_out.Linear = true;
251    fext{1} = 'df';
252    fext{2} = 'sf';
253   
254    %% Calculate T1-maps ----------------------------
255  case {'t1_ir','t1_sr','t1_3p','t1ir','t1sr','t13p'}
256   
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
264       
265    % Fit T1 maps
266    map_out = l_T1_map(datamtx,fit_vals,maptype,Dat);
267    map_out.Type = maptype;
268    map_out.Linear = false;%Dat.linear;
269        map_out.MaxIter = Dat.MaxIter;
270    fext{1} = 't1';
271    fext{2} = 's1';
272   
273    %% Calculate T1rho-maps ----------------------------
274  case {'t1r','t1rho'}
275   
276        maptype = 'T1r';
277       
278    % Fit T1r maps
279    map_out = l_T2_map(datamtx,fit_vals,maptype,Dat);
280    map_out.Type = 'T1r';
281    map_out.Linear = Dat.linear;
282    fext{1} = 't1r';
283    fext{2} = 's1r';
284   
285        %% Calculate T2rho-maps ----------------------------
286  case {'t2r','t2rho'}
287   
288        maptype = 'T2r';
289       
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';
296       
297       
298    %% Calculate T2 maps ----------------------------
299  case {'t2','r2'}
300   
301        if strcmpi(maptype,'t2')
302          % Fit T2 maps
303          map_out = l_T2_map(datamtx,fit_vals,maptype,Dat);
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
310          map_out = l_T2__map(datamtx,fit_vals,maptype,Dat);
311          map_out.Type = 'R2';
312          map_out.Linear = Dat.linear;
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
319   
320   
321    %% Calculate perfusion-maps ----------------------------
322  case {'perfusion','perf'}
323   
324   
325   
326     %% Calculate B1-maps ----------------------------
327  case 'b1'
328    maptype = 'B1';
329       
330    % Fit T1r maps
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
340    map_out.Type = maptype;
341    map_out.Linear = false;%Dat.linear;
342        map_out.MaxIter = Dat.MaxIter;
343    fext{1} = 'b1';
344    fext{2} = 'mat';
345   
346  case 'mt'
347   
348    map_out=l_MTmap(datamtx,fit_vals,Dat);
349    S0=[];
350     
351  otherwise
352    error('Unknown map type "%s"!',maptype)
353end
354
355% Add name of the parent file to the structure
356map_out.ParentFileName = ParentFileName;
357
358%% Write maps to files
359if ~isempty(Dat.filename)
360  if nargout==0
361        clear map_struct
362  end
363 
364  [fp,fn,fe]=fileparts(Dat.filename);
365  if isempty(fp)
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
377  else
378    fpath = [fp,filesep];
379  end
380  if isempty(fn)
381    fname = 'mapdata';
382  else
383    fname = fn;
384  end
385  for ii=1:size(map_out.Map,3)
386    Data = map_out.Map(:,:,ii);
387        Param = map_out;
388        Param = rmfield(Param,'Map');
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"
395        if isfield(map_out,'S0') && Dat.SaveSpinDensities
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
400  end
401else
402  map_struct = map_out;
403end
404
405
406
407%%%%%%%%%%%%%%%%%%%%%
408% Calculate T1-map
409%%%%%%%%%%%%%%%%%%%%%
410function map_out=l_T1_map(data,TI,maptype,Dat)
411% Calculate T1 relaxation times
412%
413% T1 functions:
414% Inversion recovery  :   S=abs(S0*(1-2*exp(-TI/T1)))
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
421
422%% Data size
423sz=size(data);
424
425%% Fit values
426TI = TI(:);
427
428%% Data slice index matrix
429IndMtx = reshape(1:sz(3),[length(TI) Dat.nMaps])';
430
431%% Allocate space for parameters
432map = zeros([sz(1) sz(2) Dat.nMaps]);
433S0 = zeros([sz(1) sz(2) Dat.nMaps]);
434A = zeros([sz(1) sz(2) Dat.nMaps]);
435
436% Fit options for fminsearch
437options = optimset('Display','off',...
438  'MaxIter',Dat.MaxIter);
439
440estimateInitVal = false;
441if isempty(Dat.initVal)
442  estimateInitVal = true;
443end
444
445% - Inversion recovery
446if strcmpi(maptype,'t1_ir')
447  t1_map_type = 1;
448elseif strcmpi(maptype,'t1_sr')
449  % Saturation recovery
450  t1_map_type = 2;
451else
452  % 3 Parameter fit
453  t1_map_type = 3; 
454end
455
456% Variables for waitbar
457nFits = (sz(1)*sz(2)*Dat.nMaps);
458counter = 1;
459meanTI = mean(TI);
460
461% Calculate maps
462for ii=1:Dat.nMaps
463  if Dat.wbar && ii==1
464    wbh=aedes_wbar(0,sprintf('Processing map %d/%d',ii,Dat.nMaps));
465  elseif Dat.wbar
466    aedes_wbar(counter/nFits,wbh,sprintf('Processing map %d/%d',ii,Dat.nMaps));
467  end
468 
469  for kk=1:sz(1)
470    for tt=1:sz(2)
471      if ~isempty(Dat.Mask) && ~Dat.Mask(kk,tt,ii)
472        if Dat.wbar
473          aedes_wbar(counter/nFits,wbh);
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
511        aedes_wbar(counter/nFits,wbh);
512      end
513      counter = counter+1;
514    end
515  end
516end
517if Dat.wbar
518  close(wbh)
519end
520
521map(map<0)=0;
522map(isinf(map))=0;
523map(isnan(map))=0;
524
525map_out.Map = map;
526map_out.S0 = S0;
527map_out.Angle = A;
528map_out.FitValues = TI;
529
530
531
532%%%%%%%%%%%%%%%%%%%%%%%%%%
533% Calculate B1-map
534%%%%%%%%%%%%%%%%%%%%%%%%%%
535function map_out = l_B1_map(data,fit_vals,maptype,Dat)
536
537%% Data size
538sz=size(data);
539
540%% Fit values
541fit_vals = fit_vals(:);
542
543%% Data slice index matrix
544IndMtx = reshape(1:sz(3),[length(fit_vals) Dat.nMaps])';
545
546%% Allocate space for parameters
547map = zeros([sz(1) sz(2) Dat.nMaps]);
548A = zeros([sz(1) sz(2) Dat.nMaps]);
549phase = zeros([sz(1) sz(2) Dat.nMaps]);
550
551% Fit options for fminsearch
552options = optimset('Display','off',...
553  'MaxIter',Dat.MaxIter);
554
555estimateInitVal = false;
556if isempty(Dat.initVal)
557  estimateInitVal = true;
558end
559
560
561% Variables for waitbar
562nFits = (sz(1)*sz(2)*Dat.nMaps);
563counter = 1;
564
565
566% Loop over maps
567for ii=1:Dat.nMaps
568  if Dat.wbar && ii==1
569        wbh=aedes_wbar(0,sprintf('Calculating B1-map %d/%d',ii,Dat.nMaps));
570  elseif Dat.wbar
571        aedes_wbar(counter/nFits,wbh,sprintf('Calculating map %d/%d',ii,Dat.nMaps));
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
611                fhandle= @(x) norm(data_val - x(1)*cos(x(2)*fit_vals+x(3)));
612          else
613                fhandle= @(x) norm(data_val - abs(x(1)*cos(x(2)*fit_vals+x(3))));         
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
622                aedes_wbar(counter/nFits,wbh);
623          end
624          counter = counter+1;
625        end
626  end
627end
628if Dat.wbar
629  close(wbh)
630end
631
632map(map<0)=0;
633map(isinf(map))=0;
634map(isnan(map))=0;
635
636map_out.Map = map;
637map_out.A = A;
638map_out.phase = phase;
639map_out.FitValues = fit_vals;
640
641
642%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
643% Calculate T2-, T1rho- and T2rho-maps
644%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
645function map_out=l_T2_map(data,fit_vals,maptype,Dat)
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
668done=false;
669map=[];
670S0=[];
671
672
673%% Make sure that fit values are in column vector
674fit_vals=fit_vals(:);
675
676%% Data size
677sz=size(data);
678
679%% Data slice index matrix
680IndMtx = reshape(1:sz(3),[length(fit_vals) Dat.nMaps])';
681
682map = zeros([sz(1) sz(2) Dat.nMaps]);
683S0 = zeros([sz(1) sz(2) Dat.nMaps]);
684A = zeros([sz(1) sz(2) Dat.nMaps]);
685
686% Intensity values should not be less than zero
687data(data<0)=0;
688
689%% Use linearized form (fast)
690if Dat.linear
691 
692  nFits = sz(2)*Dat.nMaps;
693  counter = 1;
694 
695  H = [ones(size(fit_vals)) -fit_vals];
696  for ii=1:Dat.nMaps
697        if Dat.wbar && ii==1
698          wbh=aedes_wbar(0,sprintf('Processing map %d/%d',ii,Dat.nMaps));
699        elseif Dat.wbar
700          aedes_wbar(counter/nFits,wbh,sprintf('Processing map %d/%d',ii,Dat.nMaps));
701        end
702       
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;
715          counter = counter+1;
716          if Dat.wbar
717                aedes_wbar(counter/nFits,wbh);
718          end
719    end
720  end
721  map(map<0)=0;
722  map(isinf(map))=0;
723  map(isnan(map))=0;
724else
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;
737  meanFV = mean(fit_vals);
738 
739  for ii=1:Dat.nMaps
740        if Dat.wbar && ii==1
741          wbh=aedes_wbar(0,sprintf('Processing map %d/%d',ii,Dat.nMaps));
742        elseif Dat.wbar
743          aedes_wbar(counter/nFits,wbh,sprintf('Processing map %d/%d',ii,Dat.nMaps));
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];
756                  init_val = [max(z); meanFV];
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);
763                fhandle = @(x) norm(z - x(1)*exp(-fit_vals./x(2)));
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
771                  aedes_wbar(counter/nFits,wbh);
772                end
773                counter = counter+1;
774          end
775        end
776  end
777 
778end
779if Dat.wbar
780  close(wbh)
781end
782
783
784map_out.Map = map;
785map_out.S0 = S0;
786map_out.FitValues = fit_vals;
787
788%%%%%%%%%%%%%%%%%%%%%%%%%%
789% Calculate Perfusion-map
790%%%%%%%%%%%%%%%%%%%%%%%%%%
791function l_Perfusion_map()
792
793
794
795%%%%%%%%%%%%%%%%%%%%%%%%%%
796% Calculate Diffusion-map
797%%%%%%%%%%%%%%%%%%%%%%%%%%
798function l_Diffusion_map()
799
800
801
802%%%%%%%%%%%%%%%%%%%%%%%%%%
803% Calculate ADC-map
804%%%%%%%%%%%%%%%%%%%%%%%%%%
805function map_out=l_ADC_map(data,b,opt,Dat)
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
817if ~exist('opt','var')
818  % Default to linearized fit
819  opt = '';
820end
821
822b=b(:);
823
824%% Data size
825sz=size(data);
826
827%% Data slice index matrix
828IndMtx = reshape(1:sz(3),[length(b) Dat.nMaps])';
829
830ADC = zeros([sz(1) sz(2) Dat.nMaps]);
831S0 = zeros([sz(1) sz(2) Dat.nMaps]);
832
833% Intensity values should not be less than zero
834data(data<0)=0;
835
836%% Use linearized form (fast)
837if isempty(opt) || strcmpi(opt,'linear')
838  H = [ones(size(b)) -b];
839  for ii=1:Dat.nMaps
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;
851else
852  %% Fit using nonlinear LS (slow, but more accurate)
853 
854end
855
856map_out.Map = ADC;
857map_out.S0 = S0;
858map_out.FitValues = b;
859
860
861
862
863
864%%%%%%%%%%%%%%%%%%%%%%%%%%
865% Calculate MT-maps
866%%%%%%%%%%%%%%%%%%%%%%%%%%
867function [map]=l_MTmap(data,fit_vals,nMaps)
868
869sz=size(data);
870lims = [-700 700];
871
872
873
874[sorted_vals,ind]=sort(fit_vals);
875[mx,mx_ind]=max(abs(sorted_vals));
876sorted_vals(mx_ind)=[];
877sorted_vals = sorted_vals(3:end-2);
878
879
880%% Allocate space for map
881map=zeros(sz(1),sz(2));
882
883lim1 = find(sorted_vals==lims(1));
884lim2 = find(sorted_vals==lims(2));
885
886
887if false
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    df=diff([vox_data(lim1) vox_data(lim2)]);
898    map(ii,kk)=df;
899  end
900end
901
902else
903
904zind=find(sorted_vals==0);
905fit=sorted_vals(zind-3:zind+3);
906fit=fit(:);
907%fit(4)=[];
908H = [fit.^2 fit ones(size(fit))];
909
910for 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))
976end
977
978end
Note: See TracBrowser for help on using the repository browser.

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