source: an2_fitmaps.m @ 44

Last change on this file since 44 was 42, checked in by tjniskan, 11 years ago
  • Rewrote the overlay feature of Aedes to work with absolute rather

than relative values.

  • The export tool can now also print overlays and ROIs with

transparencies. In these cases, however, it has to use OpenGL renderer
which is more unreliable in terms of print quality than other renderers.

  • Added a license notification that prints into the command window

every time Aedes is started. This has to do with the coming
"internationalization" of Aedes. The notification can be suppressed y
typing setpref('Aedes','ShowLicenseAtStartUp?',false) into the command
window.

  • Added some further error checking in an2_fitmaps.m
  • Fixed a minor bug in an2_write_nifti.m

M an2_export_gui.m
M an2_revision.m
M an2_fitmaps.m
M aedes.m
M an2_helpabout.m
M an2_write_nifti.m

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

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