source: an2_fitmaps.m @ 41

Last change on this file since 41 was 39, checked in by tjniskan, 11 years ago
  • Added "Recent Files" quick access menu in File-menu
  • Fixed a bug in plugins/copy_data_to_workspace.m that made writing a

new name for a variable difficult

  • Fixed a bug in the file prefix editbox in an2_export_gui.m
  • Changed references to Analyze 2.0 in all license notices to Aedes.

The name change should now be complete...

M an2_export_gui.m
M an2_cellsprintf.m
M an2_calc_wait.m
M an2_check_file_exist.m
M an2_iseven.m
M an2_cellwrite.m
M an2_wbar.m
M an2_rot3d.m
M an2_readfdf.m
M an2_revision.m
M an2_viewprocpar.m
M an2_checkcthdr.m
M an2_readprocpar.m
M an2_fitmaps.m
M an2_read_nifti.m
M an2_data_read.m
M an2_resviewer.m
M an2_maptool.m
M aedes.m
M an2_res2table.m
M an2_copy_roi.m
M plugins/save_roi_as_mask.m
M plugins/write_difference_images.m
M plugins/plot_profile.m
M plugins/calculate_t2_map.m
M plugins/calculate_t1r_map.m
M plugins/view_kspace.m
M plugins/copy_data_to_workspace.m
M plugins/take_snapshot.m
M an2_inputdlg.m
M an2_roi_copy_gui.m
M an2_readctdata.m
M an2_readfid.m
M an2_readfidprefs.m
M an2_readtab.m
M an2_check_updates.m
M an2_killfigs.m
M an2_roi_stats.m
M an2_saveres.m
M an2_rotateflip.m
M an2_juigetfiles.m
M an2_gui_defaults.m
M an2_editstack.m
M an2_update.m
M an2_write_nifti.m

File size: 22.8 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
205else
206  Dat.nMaps = sz(3)/length(fit_vals);
207 
208  %% NaN:s in fit_vals means that the corresponding slices will be omitted
209  %% from the calculations
210  indNans = find(isnan(fit_vals));
211  if not(isempty(indNans))
212        l=length(fit_vals);
213        ind = repmat(indNans,1,Dat.nMaps)+reshape(repmat(l.*[0:Dat.nMaps-1],length(indNans),1),[],1)';
214        datamtx(:,:,ind) = [];
215        fit_vals(indNans)=[];
216        sz=size(datamtx);
217  end
218 
219end
220
221
222switch lower(maptype)
223  %% Calculate ADC-maps ----------------------------
224  case {'adc','df','diffusion'}
225   
226    % Fit ADC map
227    map_out=l_ADC_map(datamtx,fit_vals,'',Dat);
228    map_out.Type = 'ADC';
229    map_out.Linear = true;
230    fext{1} = 'df';
231    fext{2} = 'sf';
232   
233    %% Calculate T1-maps ----------------------------
234  case {'t1_ir','t1_sr','t1_3p','t1ir','t1sr','t13p'}
235   
236        if any(strcmpi(maptype,{'t1_ir','t1ir'}))
237          maptype = 'T1_IR';
238        elseif any(strcmpi(maptype,{'t1_sr','t1sr'}))
239          maptype = 'T1_SR';
240        else
241          maptype = 'T1_3P';
242        end
243       
244    % Fit T1 maps
245    map_out = l_T1_map(datamtx,fit_vals,maptype,Dat);
246    map_out.Type = maptype;
247    map_out.Linear = false;%Dat.linear;
248        map_out.MaxIter = Dat.MaxIter;
249    fext{1} = 't1';
250    fext{2} = 's1';
251   
252    %% Calculate T1rho-maps ----------------------------
253  case {'t1r','t1rho'}
254   
255        maptype = 'T1r';
256       
257    % Fit T1r maps
258    map_out = l_T2_map(datamtx,fit_vals,maptype,Dat);
259    map_out.Type = 'T1r';
260    map_out.Linear = Dat.linear;
261    fext{1} = 't1r';
262    fext{2} = 's1r';
263   
264        %% Calculate T2rho-maps ----------------------------
265  case {'t2r','t2rho'}
266   
267        maptype = 'T2r';
268       
269        % Fit T2r maps
270        map_out = l_T2_map(datamtx,fit_vals,maptype,Dat);
271    map_out.Type = 'T2r';
272    map_out.Linear = Dat.linear;
273    fext{1} = 't2r';
274    fext{2} = 's2r';
275       
276       
277    %% Calculate T2 maps ----------------------------
278  case {'t2','r2'}
279   
280        if strcmpi(maptype,'t2')
281          % Fit T2 maps
282          map_out = l_T2_map(datamtx,fit_vals,maptype,Dat);
283          map_out.Type = 'T2';
284          map_out.Linear = true;
285          fext{1} = 't2';
286          fext{2} = 's2';
287        elseif strcmpi(maptype,'r2')
288          % Fit R2 maps
289          map_out = l_T2__map(datamtx,fit_vals,maptype,Dat);
290          map_out.Type = 'R2';
291          map_out.Linear = Dat.linear;
292          map_out.Map = 1./map_out.Map;
293          map_out.Map(isinf(map_out.Map))=0;
294          map_out.Map(isnan(map_out.Map))=0;
295          fext{1} = 'r2';
296          fext{2} = 's2';
297        end
298   
299   
300    %% Calculate perfusion-maps ----------------------------
301  case {'perfusion','perf'}
302   
303   
304   
305     %% Calculate B1-maps ----------------------------
306  case 'b1'
307    maptype = 'B1';
308       
309    % Fit T1r maps
310    map_out = l_B1_map(datamtx,fit_vals,maptype,Dat);
311    map_out.Type = maptype;
312    map_out.Linear = false;%Dat.linear;
313        map_out.MaxIter = Dat.MaxIter;
314    fext{1} = 'mat';
315    fext{2} = 'mat';
316   
317  case 'mt'
318   
319    map_out=l_MTmap(datamtx,fit_vals,Dat);
320    S0=[];
321     
322  otherwise
323    error('Unknown map type "%s"!',maptype)
324end
325
326% Add name of the parent file to the structure
327map_out.ParentFileName = ParentFileName;
328
329%% Write maps to files
330if ~isempty(Dat.filename)
331  if nargout==0
332        clear map_struct
333  end
334 
335  [fp,fn,fe]=fileparts(Dat.filename);
336  if isempty(fp)
337        if ~isempty(ParentFileName)
338          if iscell(ParentFileName)
339                [p,n,e]=fileparts(ParentFileName{1});
340                fpath = [p,filesep];
341          else
342                [p,n,e]=fileparts(ParentFileName);
343                fpath = [p,filesep];
344          end
345        else
346          fpath = [pwd,filesep];
347        end
348  else
349    fpath = [fp,filesep];
350  end
351  if isempty(fn)
352    fname = 'mapdata';
353  else
354    fname = fn;
355  end
356  for ii=1:size(map_out.Map,3)
357    Data = map_out.Map(:,:,ii);
358        map_out = rmfield(map_out,'Map');
359    Param = map_out;
360        %Param.Type = map_out.Type;
361    %Param.Linear = map_out.Linear;
362    %Param.FitValues = map_out.FitValues;
363        %Param.ParentFileName = map_out.ParentFileName;
364   
365    % Save map
366    filename = sprintf('%s%s_%03d.%s',fpath,fname,ii,fext{1});
367    save(filename,'Data','Param','-mat')
368   
369    % Save "spin densities"
370        if isfield(map_out,'S0') && Dat.SaveSpinDensities
371          Data = map_out.S0(:,:,ii);
372          filename = sprintf('%s%s_%03d.%s',fpath,fname,ii,fext{2});
373          save(filename,'Data','Param','-mat')
374        end
375  end
376else
377  map_struct = map_out;
378end
379
380
381
382%%%%%%%%%%%%%%%%%%%%%
383% Calculate T1-map
384%%%%%%%%%%%%%%%%%%%%%
385function map_out=l_T1_map(data,TI,maptype,Dat)
386% Calculate T1 relaxation times
387%
388% T1 functions:
389% Inversion recovery  :   S=abs(S0*(1-2*exp(-TI/T1)))
390% Saturation recovery : S=S0*(1-1*exp(-TI/T1))
391% 3 Parameter fit     : S=S0*(1-A*exp(-TI/T1))
392%
393% where:
394% S = measured signal, S0 = "spin density"
395% TI = inversion time, T1 = T1 relaxation time
396
397%% Data size
398sz=size(data);
399
400%% Fit values
401TI = TI(:);
402
403%% Data slice index matrix
404IndMtx = reshape(1:sz(3),[length(TI) Dat.nMaps])';
405
406%% Allocate space for parameters
407map = zeros([sz(1) sz(2) Dat.nMaps]);
408S0 = zeros([sz(1) sz(2) Dat.nMaps]);
409A = zeros([sz(1) sz(2) Dat.nMaps]);
410
411% Fit options for fminsearch
412options = optimset('Display','off',...
413  'MaxIter',Dat.MaxIter);
414
415estimateInitVal = false;
416if isempty(Dat.initVal)
417  estimateInitVal = true;
418end
419
420% - Inversion recovery
421if strcmpi(maptype,'t1_ir')
422  t1_map_type = 1;
423elseif strcmpi(maptype,'t1_sr')
424  % Saturation recovery
425  t1_map_type = 2;
426else
427  % 3 Parameter fit
428  t1_map_type = 3; 
429end
430
431% Variables for waitbar
432nFits = (sz(1)*sz(2)*Dat.nMaps);
433counter = 1;
434
435% Calculate maps
436for ii=1:Dat.nMaps
437  if Dat.wbar && ii==1
438        wbh=an2_wbar(0,sprintf('Processing map %d/%d',ii,Dat.nMaps));
439  elseif Dat.wbar
440        an2_wbar(counter/nFits,wbh,sprintf('Processing map %d/%d',ii,Dat.nMaps));
441  end
442
443  for kk=1:sz(1)
444        for tt=1:sz(2)
445          % Data values
446          data_val = squeeze(data(kk,tt,IndMtx(ii,:)));
447          data_val = data_val(:);
448         
449          % Initial values for the fit
450          if estimateInitVal
451                if t1_map_type == 3
452                  init_val = [mean(data_val); mean(TI); 2];
453                else
454                  init_val = [mean(data_val); mean(TI)];
455                end
456          else
457                init_val = Dat.initVal;
458          end
459         
460          % Nelder-Mead simplex iteration
461          if t1_map_type == 1
462                fhandle = @(x) sum((data_val - abs(x(1)*(1-2*exp(-TI./x(2))))).^2);
463          elseif t1_map_type == 2
464                fhandle = @(x) sum((data_val - abs(x(1)*(1-1*exp(-TI./x(2))))).^2);
465          else
466                fhandle = @(x) sum((data_val - abs(x(1)*(1-x(3)*exp(-TI./x(2))))).^2);
467          end
468          th = fminsearch(fhandle,init_val,options);
469         
470          S0(kk,tt,ii) = th(1);
471          map(kk,tt,ii) = th(2);
472          if t1_map_type == 3
473                A(kk,tt,ii) = th(3);
474          end
475         
476          if Dat.wbar
477                an2_wbar(counter/nFits,wbh);
478          end
479          counter = counter+1;
480        end
481  end
482end
483if Dat.wbar
484  close(wbh)
485end
486 
487map(map<0)=0;
488map(isinf(map))=0;
489map(isnan(map))=0;
490
491map_out.Map = map;
492map_out.S0 = S0;
493map_out.Angle = A;
494map_out.FitValues = TI;
495
496
497
498%%%%%%%%%%%%%%%%%%%%%%%%%%
499% Calculate B1-map
500%%%%%%%%%%%%%%%%%%%%%%%%%%
501function map_out = l_B1_map(data,fit_vals,maptype,Dat)
502
503%% Data size
504sz=size(data);
505
506%% Fit values
507fit_vals = fit_vals(:);
508
509%% Data slice index matrix
510IndMtx = reshape(1:sz(3),[length(fit_vals) Dat.nMaps])';
511
512%% Allocate space for parameters
513map = zeros([sz(1) sz(2) Dat.nMaps]);
514A = zeros([sz(1) sz(2) Dat.nMaps]);
515phase = zeros([sz(1) sz(2) Dat.nMaps]);
516
517% Fit options for fminsearch
518options = optimset('Display','off',...
519  'MaxIter',Dat.MaxIter);
520
521estimateInitVal = false;
522if isempty(Dat.initVal)
523  estimateInitVal = true;
524end
525
526
527% Variables for waitbar
528nFits = (sz(1)*sz(2)*Dat.nMaps);
529counter = 1;
530
531
532% Loop over maps
533for ii=1:Dat.nMaps
534  if Dat.wbar && ii==1
535        wbh=an2_wbar(0,sprintf('Calculating B1-map %d/%d',ii,Dat.nMaps));
536  elseif Dat.wbar
537        an2_wbar(counter/nFits,wbh,sprintf('Calculating map %d/%d',ii,Dat.nMaps));
538  end
539 
540  for tt=1:sz(1)
541        for kk=1:sz(2)
542          % Data values and initial values for the iteration
543          if Dat.UseComplexData
544                data_val = squeeze(data(kk,tt,IndMtx(ii,:)));
545                data_val = real(data_val);
546                data_val = data_val(:);
547               
548                if estimateInitVal
549                  val1=length(find(diff(data_val<0)<0));
550                  val2=length(find(diff(data_val<0)>0));
551                  tot_t = fit_vals(end)-fit_vals(1);
552                  omega = (((val1+val2)/2)/tot_t)*2*pi;
553                  init_val = [2*max(data_val) omega 0.1];
554                  init_val=init_val(:);
555                else
556                  init_val = Dat.initVal;
557                end
558               
559          else
560                data_val = squeeze(data(kk,tt,IndMtx(ii,:)));
561                data_val = data_val(:);
562               
563                if estimateInitVal
564                  val1=length(find(diff(data_val<(max(data_val)*0.5))<0));
565                  val2=length(find(diff(data_val<(max(data_val)*0.5))>0));
566                  tot_t = fit_vals(end)-fit_vals(1);
567                  omega = (((val1+val2)/4)/tot_t)*2*pi*0.9;
568                  init_val = [max(data_val)-min(data_val) omega 0.1];
569                  init_val=init_val(:);
570                else
571                  init_val = Dat.initVal;
572                end
573          end
574         
575          % Function handle for fminsearch
576          if Dat.UseComplexData
577                fhandle= @(x) sum((data_val - x(1)*cos(x(2)*fit_vals+x(3))).^2);
578          else
579                fhandle= @(x) sum((data_val - abs(x(1)*cos(x(2)*fit_vals+x(3)))).^2);     
580          end
581          th = fminsearch(fhandle,init_val,options);
582         
583          map(kk,tt,ii) = abs(th(2))/(2*pi);
584          A(kk,tt,ii) = th(1);
585          phase(kk,tt,ii) = th(3);
586         
587          if Dat.wbar
588                an2_wbar(counter/nFits,wbh);
589          end
590          counter = counter+1;
591        end
592  end
593end
594if Dat.wbar
595  close(wbh)
596end
597
598map(map<0)=0;
599map(isinf(map))=0;
600map(isnan(map))=0;
601
602map_out.Map = map;
603map_out.A = A;
604map_out.phase = phase;
605map_out.FitValues = fit_vals;
606
607
608%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
609% Calculate T2-, T1rho- and T2rho-maps
610%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
611function map_out=l_T2_map(data,fit_vals,maptype,Dat)
612% Calculate T1rho relaxation times
613%
614% T1rho functions:
615%
616% Linear:    log(S)=log(S0)-SL/T1rho
617% Nonlinear: S=S0*exp(-SL/T1rho)
618%
619% where S=measured signal
620%       S0 = "spin density"
621%       SL = spinlock length
622%       T1rho = T1rho relaxation time
623%
624% T2-functions:
625%
626% Linear:    log(S)=log(S0)-te/T2
627% Nonlinear: S=S0*exp(-te/T2)
628%
629% where S=measured signal
630%       S0 = "spin density"
631%       te = echo time
632%       T2 = T2 relaxation time
633
634done=false;
635map=[];
636S0=[];
637
638
639%% Make sure that fit values are in column vector
640fit_vals=fit_vals(:);
641
642%% Data size
643sz=size(data);
644
645%% Data slice index matrix
646IndMtx = reshape(1:sz(3),[length(fit_vals) Dat.nMaps])';
647
648map = zeros([sz(1) sz(2) Dat.nMaps]);
649S0 = zeros([sz(1) sz(2) Dat.nMaps]);
650A = zeros([sz(1) sz(2) Dat.nMaps]);
651if min(min(min(data)))<=0
652  data=data+abs(min(min(min(data))))+1;
653end
654
655%% Use linearized form (fast)
656if Dat.linear
657 
658  nFits = sz(2)*Dat.nMaps;
659  counter = 1;
660 
661  H = [ones(size(fit_vals)) -fit_vals];
662  for ii=1:Dat.nMaps
663        if Dat.wbar && ii==1
664          wbh=an2_wbar(0,sprintf('Processing map %d/%d',ii,Dat.nMaps));
665        elseif Dat.wbar
666          an2_wbar(counter/nFits,wbh,sprintf('Processing map %d/%d',ii,Dat.nMaps));
667        end
668       
669    for kk=1:sz(2)
670      tmp=squeeze(data(:,kk,IndMtx(ii,:))).';
671      z=log(tmp);
672      th=H\z;
673      S0(:,kk,ii)=exp(th(1,:)');
674     
675      % Try to avoid "divide by zero" warnings
676     
677      map_tmp = th(2,:)';
678      map_tmp(map_tmp==0)=eps;
679     
680      map(:,kk,ii)=1./map_tmp;
681          counter = counter+1;
682          if Dat.wbar
683                an2_wbar(counter/nFits,wbh);
684          end
685    end
686  end
687  map(map<0)=0;
688  map(isinf(map))=0;
689  map(isnan(map))=0;
690else
691  %% Fit nonlinearly using Nelder-Mead Simplex (slow, but more accurate)
692  estimateInitVal = false;
693  if isempty(Dat.initVal)
694        estimateInitVal = true;
695  end
696 
697  % Fit options
698  options = optimset('Display','off',...
699        'MaxIter',Dat.MaxIter);
700
701  nFits = sz(1)*sz(2)*Dat.nMaps;
702  counter = 1;
703 
704  for ii=1:Dat.nMaps
705        if Dat.wbar && ii==1
706          wbh=an2_wbar(0,sprintf('Processing map %d/%d',ii,Dat.nMaps));
707        elseif Dat.wbar
708          an2_wbar(counter/nFits,wbh,sprintf('Processing map %d/%d',ii,Dat.nMaps));
709        end
710       
711        for tt=1:sz(1)
712          for kk=1:sz(2)
713               
714                % Data for the fit
715                z = squeeze(data(tt,kk,IndMtx(ii,:)));
716                z=z(:);
717               
718                % Initial values for the fit
719                if estimateInitVal
720                  %init_val = [mean(z); mean(fit_vals); 1];
721                  init_val = [mean(z); mean(fit_vals)];
722                else
723                  init_val = Dat.initVal;
724                end
725               
726               
727                %fhandle = @(x) sum((z - x(1)*exp(-fit_vals./x(2))+x(3)).^2);
728                fhandle = @(x) sum((z - x(1)*exp(-fit_vals./x(2))).^2);
729                th = fminsearch(fhandle,init_val,options);
730                 
731                S0(tt,kk,ii) = th(1);
732                map(tt,kk,ii) = th(2);
733                %A(tt,kk,ii) = th(3);
734         
735                if Dat.wbar
736                  an2_wbar(counter/nFits,wbh);
737                end
738                counter = counter+1;
739          end
740        end
741  end
742 
743end
744if Dat.wbar
745  close(wbh)
746end
747
748
749map_out.Map = map;
750map_out.S0 = S0;
751map_out.FitValues = fit_vals;
752
753%%%%%%%%%%%%%%%%%%%%%%%%%%
754% Calculate Perfusion-map
755%%%%%%%%%%%%%%%%%%%%%%%%%%
756function l_Perfusion_map()
757
758
759
760%%%%%%%%%%%%%%%%%%%%%%%%%%
761% Calculate Diffusion-map
762%%%%%%%%%%%%%%%%%%%%%%%%%%
763function l_Diffusion_map()
764
765
766
767%%%%%%%%%%%%%%%%%%%%%%%%%%
768% Calculate ADC-map
769%%%%%%%%%%%%%%%%%%%%%%%%%%
770function map_out=l_ADC_map(data,b,opt,Dat)
771% Calculate Apparent Diffusion Coefficients
772% The fitted functions:
773%
774% Linear:    log(S)=log(S0)-b*ADC
775% Nonlinear: S=S0*exp(-b*ADC)
776%
777% where S=measured signal
778%       S0 = "spin density"
779%       b = diffusion weighting
780%       ADC = apparent diffusion coefficient
781
782if ~exist('opt','var')
783  % Default to linearized fit
784  opt = '';
785end
786
787b=b(:);
788
789%% Data size
790sz=size(data);
791
792%% Data slice index matrix
793IndMtx = reshape(1:sz(3),[length(b) Dat.nMaps])';
794
795ADC = zeros([sz(1) sz(2) Dat.nMaps]);
796S0 = zeros([sz(1) sz(2) Dat.nMaps]);
797if min(min(min(data)))<=0
798  data=data+abs(min(min(min(data))))+1;
799end
800
801%% Use linearized form (fast)
802if isempty(opt) || strcmpi(opt,'linear')
803  H = [ones(size(b)) -b];
804  for ii=1:Dat.nMaps
805    for kk=1:sz(2)
806      tmp=squeeze(data(:,kk,IndMtx(ii,:))).';
807      z=log(tmp);
808      th=H\z;
809      S0(:,kk,ii)=exp(th(1,:)');
810      ADC(:,kk,ii)=th(2,:)';
811    end
812  end
813  ADC(ADC<0)=0;
814  ADC(isinf(ADC))=0;
815  ADC(isnan(ADC))=0;
816else
817  %% Fit using nonlinear LS (slow, but more accurate)
818 
819end
820
821map_out.Map = ADC;
822map_out.S0 = S0;
823map_out.FitValues = b;
824
825
826
827
828
829%%%%%%%%%%%%%%%%%%%%%%%%%%
830% Calculate MT-maps
831%%%%%%%%%%%%%%%%%%%%%%%%%%
832function [map]=l_MTmap(data,fit_vals,nMaps)
833
834sz=size(data);
835lims = [-700 700];
836
837
838
839[sorted_vals,ind]=sort(fit_vals);
840[mx,mx_ind]=max(abs(sorted_vals));
841sorted_vals(mx_ind)=[];
842sorted_vals = sorted_vals(3:end-2);
843
844
845%% Allocate space for map
846map=zeros(sz(1),sz(2));
847
848lim1 = find(sorted_vals==lims(1));
849lim2 = find(sorted_vals==lims(2));
850
851
852if false
853
854for ii=1:sz(1)
855  for kk=1:sz(2)
856    vox_data = squeeze(data(ii,kk,:));
857    vox_data=vox_data(ind);
858    vox_data=vox_data./vox_data(mx_ind);
859    vox_data(mx_ind)=[];
860    vox_data = vox_data(3:end-2);
861   
862    df=diff([vox_data(lim1) vox_data(lim2)]);
863    map(ii,kk)=df;
864  end
865end
866
867else
868
869zind=find(sorted_vals==0);
870fit=sorted_vals(zind-3:zind+3);
871fit=fit(:);
872%fit(4)=[];
873H = [fit.^2 fit ones(size(fit))];
874
875for ii=1:sz(1)
876  for kk=1:sz(2)
877    vox_data = squeeze(data(ii,kk,:));
878    vox_data=vox_data(ind);
879    vox_data=vox_data./vox_data(mx_ind);
880    vox_data(mx_ind)=[];
881    vox_data = vox_data(3:end-2);
882   
883    %% Do peak picking by fitting a polynomial
884    z=vox_data(zind-3:zind+3);
885    %z(4)=[];
886    th=H\z;
887    zz=fit(1):fit(end);
888   
889    [mn,mn_ind]=min(polyval(th,zz));
890    shift_val=zz(mn_ind);
891    %map(ii,kk)=shift_val;
892    %continue
893   
894   
895    %% Do a spline interpolation to the data
896    %interp_freq = sorted_vals(1):1:sorted_vals(end);
897   
898    interp_data = pchip(sorted_vals,vox_data,...
899                        lims-shift_val);
900   
901   
902    %% Shift frequency values
903    %interp_freq = interp_freq-shift_val;
904   
905    % Get indices to limits
906    %lim1 = find(interp_freq==lims(1));
907    %lim2 = find(interp_freq==lims(2));
908   
909    try
910      df=diff([interp_data]);
911    catch
912      df=1;
913    end
914    if abs(df)<0.5
915      map(ii,kk)=df;
916    end
917   
918   
919% $$$     if ii==70 && kk==42
920% $$$       plot(fit,z,'*-',zz,polyval(th,zz),'r')
921% $$$       shift_val
922% $$$       pause
923% $$$     end
924% $$$     
925% $$$     if abs(interp_freq(mn_ind))<200
926% $$$     
927% $$$       % Shift interpolated data
928% $$$       interp_freq = interp_freq-interp_freq(mn_ind);
929% $$$       
930% $$$       % Get indices to limits
931% $$$       lim1 = find(interp_freq==lims(1));
932% $$$       lim2 = find(interp_freq==lims(2));
933% $$$       
934% $$$       df=diff([interp_data(lim1) interp_data(lim2)]);
935% $$$       if abs(df)<0.5
936% $$$         map(ii,kk)=df;
937% $$$       end
938% $$$     end
939  end
940  %disp(num2str(ii))
941end
942
943end
Note: See TracBrowser for help on using the repository browser.

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