source: an2_fitmaps.m @ 48

Last change on this file since 48 was 47, checked in by tjniskan, 11 years ago
  • Fixed a bug in an2_fitmaps which prevented saving multislice maps
  • Fixed a bug in the handling of hidden folder under Linux in an2_juigetfiles

M an2_revision.m
M an2_fitmaps.m
M an2_juigetfiles.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, 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        Param = map_out;
362        Param = rmfield(Param,'Map');
363   
364    % Save map
365    filename = sprintf('%s%s_%03d.%s',fpath,fname,ii,fext{1});
366    save(filename,'Data','Param','-mat')
367   
368    % Save "spin densities"
369        if isfield(map_out,'S0') && Dat.SaveSpinDensities
370          Data = map_out.S0(:,:,ii);
371          filename = sprintf('%s%s_%03d.%s',fpath,fname,ii,fext{2});
372          save(filename,'Data','Param','-mat')
373        end
374  end
375else
376  map_struct = map_out;
377end
378
379
380
381%%%%%%%%%%%%%%%%%%%%%
382% Calculate T1-map
383%%%%%%%%%%%%%%%%%%%%%
384function map_out=l_T1_map(data,TI,maptype,Dat)
385% Calculate T1 relaxation times
386%
387% T1 functions:
388% Inversion recovery  :   S=abs(S0*(1-2*exp(-TI/T1)))
389% Saturation recovery : S=S0*(1-1*exp(-TI/T1))
390% 3 Parameter fit     : S=S0*(1-A*exp(-TI/T1))
391%
392% where:
393% S = measured signal, S0 = "spin density"
394% TI = inversion time, T1 = T1 relaxation time
395
396%% Data size
397sz=size(data);
398
399%% Fit values
400TI = TI(:);
401
402%% Data slice index matrix
403IndMtx = reshape(1:sz(3),[length(TI) Dat.nMaps])';
404
405%% Allocate space for parameters
406map = zeros([sz(1) sz(2) Dat.nMaps]);
407S0 = zeros([sz(1) sz(2) Dat.nMaps]);
408A = zeros([sz(1) sz(2) Dat.nMaps]);
409
410% Fit options for fminsearch
411options = optimset('Display','off',...
412  'MaxIter',Dat.MaxIter);
413
414estimateInitVal = false;
415if isempty(Dat.initVal)
416  estimateInitVal = true;
417end
418
419% - Inversion recovery
420if strcmpi(maptype,'t1_ir')
421  t1_map_type = 1;
422elseif strcmpi(maptype,'t1_sr')
423  % Saturation recovery
424  t1_map_type = 2;
425else
426  % 3 Parameter fit
427  t1_map_type = 3; 
428end
429
430% Variables for waitbar
431nFits = (sz(1)*sz(2)*Dat.nMaps);
432counter = 1;
433
434% Calculate maps
435for ii=1:Dat.nMaps
436  if Dat.wbar && ii==1
437        wbh=an2_wbar(0,sprintf('Processing map %d/%d',ii,Dat.nMaps));
438  elseif Dat.wbar
439        an2_wbar(counter/nFits,wbh,sprintf('Processing map %d/%d',ii,Dat.nMaps));
440  end
441
442  for kk=1:sz(1)
443        for tt=1:sz(2)
444          % Data values
445          data_val = squeeze(data(kk,tt,IndMtx(ii,:)));
446          data_val = data_val(:);
447         
448          % Initial values for the fit
449          if estimateInitVal
450                if t1_map_type == 3
451                  init_val = [mean(data_val); mean(TI); 2];
452                else
453                  init_val = [mean(data_val); mean(TI)];
454                end
455          else
456                init_val = Dat.initVal;
457          end
458         
459          % Nelder-Mead simplex iteration
460          if t1_map_type == 1
461                fhandle = @(x) sum((data_val - abs(x(1)*(1-2*exp(-TI./x(2))))).^2);
462          elseif t1_map_type == 2
463                fhandle = @(x) sum((data_val - abs(x(1)*(1-1*exp(-TI./x(2))))).^2);
464          else
465                fhandle = @(x) sum((data_val - abs(x(1)*(1-x(3)*exp(-TI./x(2))))).^2);
466          end
467          th = fminsearch(fhandle,init_val,options);
468         
469          S0(kk,tt,ii) = th(1);
470          map(kk,tt,ii) = th(2);
471          if t1_map_type == 3
472                A(kk,tt,ii) = th(3);
473          end
474         
475          if Dat.wbar
476                an2_wbar(counter/nFits,wbh);
477          end
478          counter = counter+1;
479        end
480  end
481end
482if Dat.wbar
483  close(wbh)
484end
485 
486map(map<0)=0;
487map(isinf(map))=0;
488map(isnan(map))=0;
489
490map_out.Map = map;
491map_out.S0 = S0;
492map_out.Angle = A;
493map_out.FitValues = TI;
494
495
496
497%%%%%%%%%%%%%%%%%%%%%%%%%%
498% Calculate B1-map
499%%%%%%%%%%%%%%%%%%%%%%%%%%
500function map_out = l_B1_map(data,fit_vals,maptype,Dat)
501
502%% Data size
503sz=size(data);
504
505%% Fit values
506fit_vals = fit_vals(:);
507
508%% Data slice index matrix
509IndMtx = reshape(1:sz(3),[length(fit_vals) Dat.nMaps])';
510
511%% Allocate space for parameters
512map = zeros([sz(1) sz(2) Dat.nMaps]);
513A = zeros([sz(1) sz(2) Dat.nMaps]);
514phase = zeros([sz(1) sz(2) Dat.nMaps]);
515
516% Fit options for fminsearch
517options = optimset('Display','off',...
518  'MaxIter',Dat.MaxIter);
519
520estimateInitVal = false;
521if isempty(Dat.initVal)
522  estimateInitVal = true;
523end
524
525
526% Variables for waitbar
527nFits = (sz(1)*sz(2)*Dat.nMaps);
528counter = 1;
529
530
531% Loop over maps
532for ii=1:Dat.nMaps
533  if Dat.wbar && ii==1
534        wbh=an2_wbar(0,sprintf('Calculating B1-map %d/%d',ii,Dat.nMaps));
535  elseif Dat.wbar
536        an2_wbar(counter/nFits,wbh,sprintf('Calculating map %d/%d',ii,Dat.nMaps));
537  end
538 
539  for tt=1:sz(1)
540        for kk=1:sz(2)
541          % Data values and initial values for the iteration
542          if Dat.UseComplexData
543                data_val = squeeze(data(kk,tt,IndMtx(ii,:)));
544                data_val = real(data_val);
545                data_val = data_val(:);
546               
547                if estimateInitVal
548                  val1=length(find(diff(data_val<0)<0));
549                  val2=length(find(diff(data_val<0)>0));
550                  tot_t = fit_vals(end)-fit_vals(1);
551                  omega = (((val1+val2)/2)/tot_t)*2*pi;
552                  init_val = [2*max(data_val) omega 0.1];
553                  init_val=init_val(:);
554                else
555                  init_val = Dat.initVal;
556                end
557               
558          else
559                data_val = squeeze(data(kk,tt,IndMtx(ii,:)));
560                data_val = data_val(:);
561               
562                if estimateInitVal
563                  val1=length(find(diff(data_val<(max(data_val)*0.5))<0));
564                  val2=length(find(diff(data_val<(max(data_val)*0.5))>0));
565                  tot_t = fit_vals(end)-fit_vals(1);
566                  omega = (((val1+val2)/4)/tot_t)*2*pi*0.9;
567                  init_val = [max(data_val)-min(data_val) omega 0.1];
568                  init_val=init_val(:);
569                else
570                  init_val = Dat.initVal;
571                end
572          end
573         
574          % Function handle for fminsearch
575          if Dat.UseComplexData
576                fhandle= @(x) sum((data_val - x(1)*cos(x(2)*fit_vals+x(3))).^2);
577          else
578                fhandle= @(x) sum((data_val - abs(x(1)*cos(x(2)*fit_vals+x(3)))).^2);     
579          end
580          th = fminsearch(fhandle,init_val,options);
581         
582          map(kk,tt,ii) = abs(th(2))/(2*pi);
583          A(kk,tt,ii) = th(1);
584          phase(kk,tt,ii) = th(3);
585         
586          if Dat.wbar
587                an2_wbar(counter/nFits,wbh);
588          end
589          counter = counter+1;
590        end
591  end
592end
593if Dat.wbar
594  close(wbh)
595end
596
597map(map<0)=0;
598map(isinf(map))=0;
599map(isnan(map))=0;
600
601map_out.Map = map;
602map_out.A = A;
603map_out.phase = phase;
604map_out.FitValues = fit_vals;
605
606
607%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
608% Calculate T2-, T1rho- and T2rho-maps
609%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
610function map_out=l_T2_map(data,fit_vals,maptype,Dat)
611% Calculate T1rho relaxation times
612%
613% T1rho functions:
614%
615% Linear:    log(S)=log(S0)-SL/T1rho
616% Nonlinear: S=S0*exp(-SL/T1rho)
617%
618% where S=measured signal
619%       S0 = "spin density"
620%       SL = spinlock length
621%       T1rho = T1rho relaxation time
622%
623% T2-functions:
624%
625% Linear:    log(S)=log(S0)-te/T2
626% Nonlinear: S=S0*exp(-te/T2)
627%
628% where S=measured signal
629%       S0 = "spin density"
630%       te = echo time
631%       T2 = T2 relaxation time
632
633done=false;
634map=[];
635S0=[];
636
637
638%% Make sure that fit values are in column vector
639fit_vals=fit_vals(:);
640
641%% Data size
642sz=size(data);
643
644%% Data slice index matrix
645IndMtx = reshape(1:sz(3),[length(fit_vals) Dat.nMaps])';
646
647map = zeros([sz(1) sz(2) Dat.nMaps]);
648S0 = zeros([sz(1) sz(2) Dat.nMaps]);
649A = zeros([sz(1) sz(2) Dat.nMaps]);
650if min(min(min(data)))<=0
651  data=data+abs(min(min(min(data))))+1;
652end
653
654%% Use linearized form (fast)
655if Dat.linear
656 
657  nFits = sz(2)*Dat.nMaps;
658  counter = 1;
659 
660  H = [ones(size(fit_vals)) -fit_vals];
661  for ii=1:Dat.nMaps
662        if Dat.wbar && ii==1
663          wbh=an2_wbar(0,sprintf('Processing map %d/%d',ii,Dat.nMaps));
664        elseif Dat.wbar
665          an2_wbar(counter/nFits,wbh,sprintf('Processing map %d/%d',ii,Dat.nMaps));
666        end
667       
668    for kk=1:sz(2)
669      tmp=squeeze(data(:,kk,IndMtx(ii,:))).';
670      z=log(tmp);
671      th=H\z;
672      S0(:,kk,ii)=exp(th(1,:)');
673     
674      % Try to avoid "divide by zero" warnings
675     
676      map_tmp = th(2,:)';
677      map_tmp(map_tmp==0)=eps;
678     
679      map(:,kk,ii)=1./map_tmp;
680          counter = counter+1;
681          if Dat.wbar
682                an2_wbar(counter/nFits,wbh);
683          end
684    end
685  end
686  map(map<0)=0;
687  map(isinf(map))=0;
688  map(isnan(map))=0;
689else
690  %% Fit nonlinearly using Nelder-Mead Simplex (slow, but more accurate)
691  estimateInitVal = false;
692  if isempty(Dat.initVal)
693        estimateInitVal = true;
694  end
695 
696  % Fit options
697  options = optimset('Display','off',...
698        'MaxIter',Dat.MaxIter);
699
700  nFits = sz(1)*sz(2)*Dat.nMaps;
701  counter = 1;
702 
703  for ii=1:Dat.nMaps
704        if Dat.wbar && ii==1
705          wbh=an2_wbar(0,sprintf('Processing map %d/%d',ii,Dat.nMaps));
706        elseif Dat.wbar
707          an2_wbar(counter/nFits,wbh,sprintf('Processing map %d/%d',ii,Dat.nMaps));
708        end
709       
710        for tt=1:sz(1)
711          for kk=1:sz(2)
712               
713                % Data for the fit
714                z = squeeze(data(tt,kk,IndMtx(ii,:)));
715                z=z(:);
716               
717                % Initial values for the fit
718                if estimateInitVal
719                  %init_val = [mean(z); mean(fit_vals); 1];
720                  init_val = [mean(z); mean(fit_vals)];
721                else
722                  init_val = Dat.initVal;
723                end
724               
725               
726                %fhandle = @(x) sum((z - x(1)*exp(-fit_vals./x(2))+x(3)).^2);
727                fhandle = @(x) sum((z - x(1)*exp(-fit_vals./x(2))).^2);
728                th = fminsearch(fhandle,init_val,options);
729                 
730                S0(tt,kk,ii) = th(1);
731                map(tt,kk,ii) = th(2);
732                %A(tt,kk,ii) = th(3);
733         
734                if Dat.wbar
735                  an2_wbar(counter/nFits,wbh);
736                end
737                counter = counter+1;
738          end
739        end
740  end
741 
742end
743if Dat.wbar
744  close(wbh)
745end
746
747
748map_out.Map = map;
749map_out.S0 = S0;
750map_out.FitValues = fit_vals;
751
752%%%%%%%%%%%%%%%%%%%%%%%%%%
753% Calculate Perfusion-map
754%%%%%%%%%%%%%%%%%%%%%%%%%%
755function l_Perfusion_map()
756
757
758
759%%%%%%%%%%%%%%%%%%%%%%%%%%
760% Calculate Diffusion-map
761%%%%%%%%%%%%%%%%%%%%%%%%%%
762function l_Diffusion_map()
763
764
765
766%%%%%%%%%%%%%%%%%%%%%%%%%%
767% Calculate ADC-map
768%%%%%%%%%%%%%%%%%%%%%%%%%%
769function map_out=l_ADC_map(data,b,opt,Dat)
770% Calculate Apparent Diffusion Coefficients
771% The fitted functions:
772%
773% Linear:    log(S)=log(S0)-b*ADC
774% Nonlinear: S=S0*exp(-b*ADC)
775%
776% where S=measured signal
777%       S0 = "spin density"
778%       b = diffusion weighting
779%       ADC = apparent diffusion coefficient
780
781if ~exist('opt','var')
782  % Default to linearized fit
783  opt = '';
784end
785
786b=b(:);
787
788%% Data size
789sz=size(data);
790
791%% Data slice index matrix
792IndMtx = reshape(1:sz(3),[length(b) Dat.nMaps])';
793
794ADC = zeros([sz(1) sz(2) Dat.nMaps]);
795S0 = zeros([sz(1) sz(2) Dat.nMaps]);
796if min(min(min(data)))<=0
797  data=data+abs(min(min(min(data))))+1;
798end
799
800%% Use linearized form (fast)
801if isempty(opt) || strcmpi(opt,'linear')
802  H = [ones(size(b)) -b];
803  for ii=1:Dat.nMaps
804    for kk=1:sz(2)
805      tmp=squeeze(data(:,kk,IndMtx(ii,:))).';
806      z=log(tmp);
807      th=H\z;
808      S0(:,kk,ii)=exp(th(1,:)');
809      ADC(:,kk,ii)=th(2,:)';
810    end
811  end
812  ADC(ADC<0)=0;
813  ADC(isinf(ADC))=0;
814  ADC(isnan(ADC))=0;
815else
816  %% Fit using nonlinear LS (slow, but more accurate)
817 
818end
819
820map_out.Map = ADC;
821map_out.S0 = S0;
822map_out.FitValues = b;
823
824
825
826
827
828%%%%%%%%%%%%%%%%%%%%%%%%%%
829% Calculate MT-maps
830%%%%%%%%%%%%%%%%%%%%%%%%%%
831function [map]=l_MTmap(data,fit_vals,nMaps)
832
833sz=size(data);
834lims = [-700 700];
835
836
837
838[sorted_vals,ind]=sort(fit_vals);
839[mx,mx_ind]=max(abs(sorted_vals));
840sorted_vals(mx_ind)=[];
841sorted_vals = sorted_vals(3:end-2);
842
843
844%% Allocate space for map
845map=zeros(sz(1),sz(2));
846
847lim1 = find(sorted_vals==lims(1));
848lim2 = find(sorted_vals==lims(2));
849
850
851if false
852
853for ii=1:sz(1)
854  for kk=1:sz(2)
855    vox_data = squeeze(data(ii,kk,:));
856    vox_data=vox_data(ind);
857    vox_data=vox_data./vox_data(mx_ind);
858    vox_data(mx_ind)=[];
859    vox_data = vox_data(3:end-2);
860   
861    df=diff([vox_data(lim1) vox_data(lim2)]);
862    map(ii,kk)=df;
863  end
864end
865
866else
867
868zind=find(sorted_vals==0);
869fit=sorted_vals(zind-3:zind+3);
870fit=fit(:);
871%fit(4)=[];
872H = [fit.^2 fit ones(size(fit))];
873
874for ii=1:sz(1)
875  for kk=1:sz(2)
876    vox_data = squeeze(data(ii,kk,:));
877    vox_data=vox_data(ind);
878    vox_data=vox_data./vox_data(mx_ind);
879    vox_data(mx_ind)=[];
880    vox_data = vox_data(3:end-2);
881   
882    %% Do peak picking by fitting a polynomial
883    z=vox_data(zind-3:zind+3);
884    %z(4)=[];
885    th=H\z;
886    zz=fit(1):fit(end);
887   
888    [mn,mn_ind]=min(polyval(th,zz));
889    shift_val=zz(mn_ind);
890    %map(ii,kk)=shift_val;
891    %continue
892   
893   
894    %% Do a spline interpolation to the data
895    %interp_freq = sorted_vals(1):1:sorted_vals(end);
896   
897    interp_data = pchip(sorted_vals,vox_data,...
898                        lims-shift_val);
899   
900   
901    %% Shift frequency values
902    %interp_freq = interp_freq-shift_val;
903   
904    % Get indices to limits
905    %lim1 = find(interp_freq==lims(1));
906    %lim2 = find(interp_freq==lims(2));
907   
908    try
909      df=diff([interp_data]);
910    catch
911      df=1;
912    end
913    if abs(df)<0.5
914      map(ii,kk)=df;
915    end
916   
917   
918% $$$     if ii==70 && kk==42
919% $$$       plot(fit,z,'*-',zz,polyval(th,zz),'r')
920% $$$       shift_val
921% $$$       pause
922% $$$     end
923% $$$     
924% $$$     if abs(interp_freq(mn_ind))<200
925% $$$     
926% $$$       % Shift interpolated data
927% $$$       interp_freq = interp_freq-interp_freq(mn_ind);
928% $$$       
929% $$$       % Get indices to limits
930% $$$       lim1 = find(interp_freq==lims(1));
931% $$$       lim2 = find(interp_freq==lims(2));
932% $$$       
933% $$$       df=diff([interp_data(lim1) interp_data(lim2)]);
934% $$$       if abs(df)<0.5
935% $$$         map(ii,kk)=df;
936% $$$       end
937% $$$     end
938  end
939  %disp(num2str(ii))
940end
941
942end
Note: See TracBrowser for help on using the repository browser.

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