source: plugins/calc_asl_cbf.m

Last change on this file was 185, checked in by tjniskan, 8 years ago
  • Added an option for calculating CBF time series from CASL data in ASL-CBF plugin.

M plugins/calc_asl_cbf.m
M aedes_revision.m

File size: 5.7 KB
Line 
1function calc_asl_cbf(DATA,ROI,AddInfo)
2% CALC_ASL_CBF - Calculate CBF maps from ASL data using T1-map data if
3% possible (Aedes plugin)
4%   
5%
6% Synopsis:
7%
8% Description:
9%
10% Examples:
11%
12% See also:
13%
14
15% This function is a part of Aedes - A graphical tool for analyzing
16% medical images
17%
18% Copyright (C) 2006 Juha-Pekka Niskanen <Juha-Pekka.Niskanen@uku.fi>
19%
20% Department of Physics, Department of Neurobiology
21% University of Kuopio, FINLAND
22%
23% This program may be used under the terms of the GNU General Public
24% License version 2.0 as published by the Free Software Foundation
25% and appearing in the file LICENSE.TXT included in the packaging of
26% this program.
27%
28% This program is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
29% WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
30
31% Defaults
32lambda = 0.9;
33t1map = 1.2;
34
35if AddInfo.isDataMixed
36  % Check that the number of images is OK
37  if rem(length(DATA),2)~=0 || length(DATA)==1
38        h=errordlg('The number of input images has to be divisible by 2!',...
39          'Error','modal');
40        return
41  end
42  controlIm = zeros([size(DATA{1}.FTDATA) length(DATA)/2]);
43  labelIm = zeros([size(DATA{1}.FTDATA) length(DATA)/2]);
44  for ii=1:2:size(controlIm,3)
45        controlIm(:,:,ii) = DATA{ii}.FTDATA;
46        labelIm(:,:,ii) = DATA{ii+1}.FTDATA;
47  end
48  for ii=1:length(DATA)
49        ParentFileName{ii} = fullfile(DATA{ii}.HDR.fpath,DATA{ii}.HDR.fname);
50  end
51else
52       
53        if strcmpi(DATA{1}.DataFormat,'bruker_reco')
54                % Assume that Bruker data is 4D
55                if rem(size(DATA{1}.FTDATA,4),2)~=0 || size(DATA{1}.FTDATA,4)==1
56                        h=errordlg('The number of input images has to be divisible by 2!',...
57                                'Error','modal');
58                        return
59                end
60               
61                controlIm = squeeze(DATA{1}.FTDATA(:,:,2,:));
62                labelIm = squeeze(DATA{1}.FTDATA(:,:,1,:));
63                ParentFileName = fullfile(DATA{1}.HDR.fpath,DATA{1}.HDR.fname);
64        else
65                % Check that the number of images is OK
66                if rem(size(DATA{1}.FTDATA,3),2)~=0 || size(DATA{1}.FTDATA,3)==1
67                        h=errordlg('The number of input images has to be divisible by 2!',...
68                                'Error','modal');
69                        return
70                end
71                if isfield(DATA{1}.PROCPAR,'ltype') && ...
72                                length(DATA{1}.PROCPAR.ltype)==size(DATA{1}.FTDATA,3)
73                        control_ind = find(strcmpi(DATA{1}.PROCPAR.ltype,'c'));
74                        label_ind = find(strcmpi(DATA{1}.PROCPAR.ltype,'l'));
75                        controlIm = DATA{1}.FTDATA(:,:,control_ind);
76                        labelIm = DATA{1}.FTDATA(:,:,label_ind);
77                        ParentFileName = fullfile(DATA{1}.HDR.fpath,DATA{1}.HDR.fname);
78                else
79                        controlIm = DATA{1}.FTDATA(:,:,1:2:size(DATA{1}.FTDATA,3));
80                        labelIm = DATA{1}.FTDATA(:,:,2:2:size(DATA{1}.FTDATA,3));
81                        ParentFileName = fullfile(DATA{1}.HDR.fpath,DATA{1}.HDR.fname);
82                end
83        end
84end
85
86% Ask for lambda
87resp = aedes_inputdlg('Lambda coefficient?','Lambda coefficient?',...
88  num2str(lambda));
89if isempty(resp)
90  % Canceled
91  return
92end
93lambda = str2num(resp{1});
94
95% Ask for T1-map
96default_path = DATA{1}.HDR.fpath;
97if isempty(default_path)
98  default_path = [pwd,filesep];
99end
100
101resp = questdlg(['Select if you want to use a constant value for T1 ',...
102  'or an existing T1-map.'],...
103  'Use constant T1 value or T1-map?',...
104  'Use T1-map','Use constant T1','Cancel',...
105  'Use T1-map');
106if isempty(resp) || strcmpi(resp,'Cancel')
107  % Canceled
108  return
109elseif strcmpi(resp,'Use constant T1')
110  % Ask for T1 value
111  resp = aedes_inputdlg('Input T1 value','Input T1 value',...
112        num2str(t1map));
113  if isempty(resp)
114        % Canceled
115        return
116  end
117  fname = '';
118  fpath = '';
119  t1map = str2num(resp{1});
120  t1val = t1map;
121else
122  [fname,fpath,findex] = uigetfile({'*.t1','T1-Files (*.t1)';...
123        '*.*','All Files (*.*)'},'Select T1-map to open',...
124        default_path);
125  if isequal(fname,0)
126        % Canceled
127        return
128  end
129  try
130        tmp=load(fullfile(fpath,fname),'-mat');
131        t1map = tmp.Data;
132  catch
133        h=errordlg({'Could not open T1-map from file',...
134          fullfile(fpath,fname)},'Error','modal');
135        return
136  end
137 
138
139  % Check that the T1-map size matches with ASL. If it doesn't, resize it to
140  % match...
141  t1map = t1map(:,:,1);
142  if ~isequal(size(t1map),[size(controlIm,1),size(controlIm,2)])
143        t1map = imresize(t1map,[size(controlIm,1),size(controlIm,2)]);
144  end
145  t1val = [];
146end
147
148% Ask if the user wants to calculate a single asl map or time series
149resp = questdlg(['Calculate an average CBF map (recommended) ',...
150        'or a CBF time series (i.e. one CBF map / image pair)?'],...
151  'Calculate an average CBF map?',...
152  'Average CBF','CBF map / image pair','Cancel',...
153  'Average CBF');
154if isempty(resp) || strcmpi(resp,'Cancel')
155        % Canceled
156        return
157end
158if strcmpi(resp,'Average CBF')
159        CBFts = false;
160else
161        CBFts = true;
162end
163
164% Convert milliseconds to seconds
165%t1map=t1map./1000;
166
167if CBFts
168        % Calculate differences and CBF maps
169        asl_map = zeros(size(controlIm));
170        cbf_map = zeros(size(controlIm));
171        for ii=1:size(controlIm,3)
172                asl_map(:,:,ii) = controlIm(:,:,ii)-labelIm(:,:,ii);
173                cbf_map(:,:,ii) = (100.*60.*lambda.*asl_map(:,:,ii))./(t1map.*2.*controlIm(:,:,ii));
174        end
175
176else
177        % Calculate difference fotr ASL
178        mean_control = mean(controlIm,3);
179        mean_label = mean(labelIm,3);
180        asl_map = mean_control-mean_label;
181       
182        % Calculate CBF and save it into a MAT-file
183        cbf_map = (100.*60.*lambda.*asl_map)./(t1map.*2.*mean_control);
184
185end
186cbf_map(isinf(cbf_map))=0;
187cbf_map(isnan(cbf_map))=0;
188
189% Ask where to save the resulting CBF map
190[fn,fp,fi] = uiputfile({'*.mat','Matlab MAT-Files (*.mat)';...
191        '*.*','All Files (*.*)'},'Save CBF-map as...',...
192        [default_path,'asl_cbf_map.mat']);
193if isequal(fname,0)
194        % Canceled
195        return
196end
197
198% Construct the map structure
199Data = cbf_map;
200Param.ParentFileName = ParentFileName;
201Param.T1FileName = fullfile(fpath,fname);
202Param.T1value = t1val;
203Param.Lambda = lambda;
204
205% Save the map
206try
207  save(fullfile(fp,fn),'Data','Param','-mat');
208catch
209  h=errordlg({'Could not save CBF-map to',...
210        fullfile(fp,fn)},'Error','modal');
211  return
212end
213
214% Open in a new Aedes window
215aedes(fullfile(fp,fn))
216
217
218
219
220
221
222
223
224
225
Note: See TracBrowser for help on using the repository browser.

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