source: plugins/fat_analysis.m @ 161

Last change on this file since 161 was 155, checked in by tjniskan, 8 years ago
  • Added support for multiple monitor systems

M aedes_inputdlg.m
M aedes_roi_copy_gui.m
A aedes_dialoglocation.m
M aedes_readfidprefs.m
M aedes_rotateflip.m
M aedes_resviewer.m
M aedes.m
M aedes_helpabout.m
M plugins/fat_analysis.m
M plugins/copy_data_to_workspace.m
M aedes_export_gui.m
M aedes_headerbrowser.m
M aedes_revision.m
M aedes_juigetfiles.m
M aedes_editstack.m

File size: 10.7 KB
Line 
1function fat_analysis(DATA,ROI,AddInfo)
2% FAT_ANALYSIS - Calculate fat percentage (Aedes plugin)
3%   
4%
5% Synopsis:
6%
7% Description:
8%
9% Examples:
10%
11% See also:
12%
13
14% This function is a part of Aedes - A graphical tool for analyzing
15% medical images
16%
17% Copyright (C) 2006 Juha-Pekka Niskanen <Juha-Pekka.Niskanen@uku.fi>
18%
19% Department of Physics, Department of Neurobiology
20% University of Kuopio, FINLAND
21%
22% This program may be used under the terms of the GNU General Public
23% License version 2.0 as published by the Free Software Foundation
24% and appearing in the file LICENSE.TXT included in the packaging of
25% this program.
26%
27% This program is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
28% WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
29
30
31%% Check if the data to be analyzed is valid
32if AddInfo.isDataMixed
33  errordlg('Sorry. This plugin does not work with mixed data.','Error',...
34    'modal');
35  return
36end
37
38%% Check if the size of the data is correct
39data_sz = size(DATA{1}.FTDATA);
40if length(data_sz)<4 || data_sz(4)~=2
41  errordlg(['The data should be a 4D-matrix containing fat and water ',...
42    'volumes in the 4th dimesion.'],'Error',...
43    'modal');
44  return
45end
46
47
48%% Warn if data does not appear to be correct but give a choice to continue
49%% anyway...
50if ~isfield(DATA{1},'PROCPAR') || ...
51    ~isfield(DATA{1}.PROCPAR,'seqfil') || ...
52    ~strcmp(DATA{1}.PROCPAR.seqfil,'ge3d_csi2')
53   
54  resp=questdlg(['Data does not seem to be measured with a valid sequence. ',...
55    'Do you want to continue anyway?'],...
56    'Invalid sequence detected',...
57    'Yes','No','No');
58  if isempty(resp) || strcmpi(resp,'no')
59    % Cancelled
60    return
61  end
62end
63
64%% Prompt for the slices that correspond to right anatomical area and have
65%% tolerable level of noise/artefacts
66canceled=false;
67invalid_input = true;
68while invalid_input
69  resp = aedes_inputdlg(['Input the appropriate slice range'],...
70    'Input slice range','1:end');
71  if isempty(resp)
72    % Cancelled
73    canceled=true;
74    invalid_input = false;
75    continue
76  end
77  resp = strrep(resp{1},'end',num2str(data_sz(3)));
78  z_range = str2num(resp);
79  if ~isempty(z_range) && isnumeric(z_range) && isreal(z_range)
80    invalid_input = false;
81  else
82    h=errordlg('Invalid range input','Invalid range input','modal');
83    uiwait(h);
84  end
85end
86if canceled
87  return
88end
89
90%% Get FOV
91if ~isfield(DATA{1},'PROCPAR') || ~isfield(DATA{1}.PROCPAR,'lro')
92  %% Prompt for FOV
93  resp = aedes_inputdlg('Input FOV in millimeters',...
94    'Input FOV','[40 40 80]');
95  if isempty(resp)
96    % Canceled
97    return
98  end
99  FOV = str2num(resp{1});
100  if isempty(FOV) || length(FOV)~=3
101    errordlg('Invalid FOV! Aborting...','Invalid FOV','modal');
102    return
103  end
104else
105  % FOV in millimeters
106  FOV = [DATA{1}.PROCPAR.lro,DATA{1}.PROCPAR.lpe,...
107    DATA{1}.PROCPAR.lpe2]*10;
108end
109
110% A brute force -method for correct DC artifacts, may not always work
111data_water = DATA{1}.FTDATA(:,:,:,2);
112data_fat = DATA{1}.FTDATA(:,:,:,1);
113[mx,DC_Ind_fat]=max(max(max(data_fat,[],2),[],3));
114[mx,DC_Ind_water]=max(max(max(data_water,[],2),[],3));
115
116% Check that the found noise peak is near the center
117if DC_Ind_fat < data_sz(1)/2-5 || DC_Ind_fat > data_sz(1)/2+5
118  warning('Largest noise peak not near the center line in fat data! Check images and results.');
119  DC_Ind_fat=data_sz(1)/2+1;
120end
121if DC_Ind_water < data_sz(1)/2-5 || DC_Ind_water > data_sz(1)/2+5
122  warning('Largest noise peak not near the center line in water data! Check images and results.');
123  DC_Ind_water=data_sz(1)/2+1;
124end
125data_fat(DC_Ind_fat,:,:) = mean(data_fat([DC_Ind_fat-1 DC_Ind_fat+1],:,:));
126data_water(DC_Ind_water,:,:) = mean(data_water([DC_Ind_water-1 DC_Ind_water+1],:,:));
127
128
129%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
130%% FAT analysis -----------------------------------------
131
132
133%% Set the default noise ROI
134% noise_fat = 2*squeeze(mean(mean(DATA{1}.FTDATA(data_sz(1)-10:data_sz(1),1:data_sz(2),:,1),1),2));
135% noise_water = 2*squeeze(mean(mean(DATA{1}.FTDATA(data_sz(1)-10:data_sz(1),1:data_sz(2),:,2),1),2));
136noise_fat = 2*squeeze(mean(mean(data_fat(1:30,1:40,:),1),2));
137noise_water = 2*squeeze(mean(mean(data_water(1:30,1:40,:),1),2));
138
139mean_fat = squeeze(mean(mean(data_fat(:,:,:),1),2));
140mean_water = squeeze(mean(mean(data_water(:,:,:),1),2));
141
142fat_th = noise_fat+mean_fat;
143water_th = noise_water+mean_water;
144fat_th = fat_th(:);
145water_th = water_th(:);
146
147% Construct matrices for threshold comparison
148fat_ind = repmat(permute(fat_th(z_range),[3 2 1]),...
149  [data_sz(1) data_sz(2) 1]);
150water_ind = repmat(permute(water_th(z_range),[3 2 1]),...
151  [data_sz(1) data_sz(2) 1]);
152
153% Calculate number of water and fat voxels
154nVoxFat = numel(find(data_fat(:,:,z_range)>fat_ind));
155nVoxWater = numel(find(data_water(:,:,z_range)>water_ind));
156
157% Calculate fat and water concentrations
158density_fat = 0.9; %g/ml;
159voxel_weight_fat = (FOV(1)/data_sz(1)*FOV(2)/data_sz(2)*FOV(3)/data_sz(3))*density_fat/1000;
160fat_weight = nVoxFat*voxel_weight_fat;
161
162density_water = 1; %g/ml;
163voxel_weight_water = (FOV(1)/data_sz(1)*FOV(2)/data_sz(2)*FOV(3)/data_sz(3))*density_water/1000;
164water_weight = nVoxWater*voxel_weight_water;
165
166% Fat-%
167fat_percent = fat_weight/(fat_weight+water_weight)*100;
168
169% Print results to the command window
170fprintf(1,'\n####################################################\n');
171fprintf(1,'FAT ANALYSIS RESULTS\n\n');
172fprintf(1,'File: %s\n',[DATA{1}.HDR.fpath,DATA{1}.HDR.fname]);
173fprintf(1,'FOV: %dx%dx%d mm\n',FOV);
174fprintf(1,'Number of fat voxels:\t%d\n',nVoxFat);
175fprintf(1,'Number of H2O voxels:\t%d\n\n',nVoxWater);
176fprintf(1,'==========================================\n');
177fprintf(1,'%s\n',fliplr(sprintf('%15s%15s%15s',...
178  fliplr('Fat prc (%)'),fliplr('H2O (mg/ml)'),fliplr('Fat (mg/ml)'))));
179fprintf(1,'------------------------------------------\n')
180fprintf(1,'%s\n',fliplr(sprintf('%15s%15s%15s',...
181  fliplr(sprintf('%.4f',fat_percent)),...
182  fliplr(sprintf('%.4f',water_weight)),fliplr(sprintf('%.4f',fat_weight)))));
183fprintf(1,'==========================================\n');
184fprintf(1,'####################################################\n');
185
186
187% Show results in uitable
188fig_h = 150;
189fig_w = 350;
190fig_location = aedes_dialoglocation([fig_w,fig_h]);
191fig_pos = [fig_location(1) fig_location(2) fig_w fig_h];
192fh = figure('units','pixels',...
193  'position',fig_pos,...
194  'Name','Fat Analysis Results',...
195  'IntegerHandle','off',...
196  'Numbertitle','off',...
197  'MenuBar','none');
198
199matlab_ver = version;
200if str2num(matlab_ver(1:3))<=7.5
201  h_uitable = uitable('parent',fh,...
202    'position',[5 5 340 55],...
203    'Data',[fat_weight water_weight fat_percent],...
204    'ColumnNames',{'Fat (mg/ml)','H2O (mg/ml)','Fat prc (%)'},...
205    'ColumnWidth',100);
206else
207  h_uitable = uitable('parent',fh,...
208    'position',[5 5 340 55],...
209    'Data',[fat_weight water_weight fat_percent],...
210    'ColumnName',{'Fat (mg/ml)','H2O (mg/ml)','Fat prc (%)'},...
211    'ColumnWidth',{100 100 100});
212end
213
214pos = get(h_uitable,'position');
215h_text = uicontrol('parent',fh,...
216  'style','text',...
217  'position',[pos(1) pos(2)+pos(4)+2 pos(3) fig_h-pos(2)-pos(4)-4],...
218  'string',...
219  {['File: ',DATA{1}.HDR.fpath,DATA{1}.HDR.fname],...
220  sprintf('FOV: %dx%dx%d mm',FOV),...
221  sprintf('Number of fat voxels: %d',nVoxFat),...
222  sprintf('Number of H2O voxels: %d',nVoxWater)},...
223  'horizontalalign','left');
224
225if ispc
226  set(h_text,'fontsize',8);
227else
228  set(h_text,'fontsize',8);
229end
230
231
232% Prompt for saving image
233resp = questdlg(['Do you also want to save images of the fat/water distributions results?'],...
234  'Save images?','Yes','No','Yes');
235if isempty(resp) || strcmpi(resp,'No')
236  % Cancelled
237  return
238end
239
240% Prompt for output directory and file name
241if ~isfield(DATA{1}.HDR,'fpath') || isempty(DATA{1}.HDR.fpath)
242  default_name = ['fat_analysis_',datestr(now,30)];
243else
244  if strcmpi(DATA{1}.HDR.fpath(end),filesep)
245    [fp,fn,fext]=fileparts(DATA{1}.HDR.fpath(1:end-1));
246  else
247    [fp,fn,fext]=fileparts(DATA{1}.HDR.fpath);
248  end
249  default_name = [fn,'_fat_analysis'];
250end
251
252try
253  default_path = getpref('FatAnalysis','ImagePath');
254catch
255  default_path = [pwd,filesep];
256end
257
258[fname,fpath,findex] = uiputfile(...
259  {'*.png','PNG Image Files (*.png)';...
260  '*.tif','TIF Image Files (*.tif)'},...
261  'Save Images to...',...
262  [default_path,default_name]);
263if isequal(fname,0) || isequal(fpath,0)
264  % Cancelled
265  return
266end
267setpref('FatAnalysis','ImagePath',fpath);
268
269% Get sw from procpar or prompt for it if it is missing
270if ~isfield(DATA{1},'PROCPAR') || ~isfield(DATA{1}.PROCPAR,'sw')
271  resp = aedes_inputdlg('Please give sw value',...
272    'Input sw value','');
273  if isempty(resp)
274    % Cancelled
275    return
276  end
277  sw = str2num(resp{1});
278  if isempty(sw) || ~isreal(sw)
279    errordlg('Invalid value for sw! Aborting...','Invalid sw value','modal');
280    return
281  end
282else
283  sw = DATA{1}.PROCPAR.sw;
284end
285
286
287% A brute force -method for correct DC artifacts, may not always work
288% data_water = DATA{1}.FTDATA(:,:,:,2);
289% data_fat_shifted = DATA{1}.FTDATA(:,:,:,1);
290% [mx,DC_Ind_fat]=max(max(DATA{1}.FTDATA(:,:,size(DATA{1}.FTDATA,3)/2,1),[],2));
291% [mx,DC_Ind_water]=max(max(DATA{1}.FTDATA(:,:,size(DATA{1}.FTDATA,3)/2,2),[],2));
292% data_fat_shifted(DC_Ind_fat,:,:) = mean(data_fat_shifted([DC_Ind_fat-1 DC_Ind_fat+1],:,:));
293% data_water(DC_Ind_water,:,:) = mean(data_water([DC_Ind_water-1 DC_Ind_water+1],:,:));
294
295%% Calculate Fat shift ------------------
296
297% Show waitbar
298[h,txt]=aedes_calc_wait({'Correcting fat shift.','This may take some time...'});
299fat_shift = 680/sw;
300
301data_fat = interp1(1:size(data_fat,1),...
302  data_fat,...
303  (1:size(data_fat,1))+fat_shift,'spline',0);
304
305% If there is a huge DC artifact in the data the spile interpolation tends
306% to go crazy near the artifacts and produce huge negative values...
307data_fat(data_fat<0) = 0;
308
309% Close waitbar
310close(h)
311% ---------------------------------------
312
313% Try to set the intensity levels to reasonable values
314%tmp_water=data_water;
315tmp_water = data_water((size(data_water,1)/2-30):(size(data_water,1)/2+30),...
316  (size(data_water,2)/2-30):(size(data_water,2)/2+30),:);
317tmp_fat = data_fat((size(data_fat,1)/2-30):(size(data_fat,1)/2+30),...
318  (size(data_fat,2)/2-30):(size(data_fat,2)/2+30),:);
319data_fat = data_fat./max(tmp_fat(:));
320data_water = data_water./max(tmp_water(:));
321data_fat(data_fat>1) = 1;
322data_water(data_water>1) = 1;
323
324% Construct an RGB image of water and fat content
325RGB_im = cat(4,data_fat,...
326  zeros(size(data_fat)),...
327  data_water);
328
329% Write png or tif images
330nImages = size(RGB_im,3);
331[fp,fn,fext]=fileparts(fname);
332if ~isempty(findex) && findex==2
333  imformat = 'TIF';
334  f_ext = '.tif';
335else
336  imformat = 'PNG';
337  f_ext = '.png';
338end
339 
340h = aedes_wbar(0,'');
341for ii=1:nImages
342  filename = sprintf([fn,'_%03d',f_ext],ii);
343  aedes_wbar(ii/nImages,h,{sprintf('Writing image %d/%d to',ii,nImages),...
344    [fpath,filename]});
345  imwrite(squeeze(RGB_im(:,:,ii,:)),[fpath,filename],imformat);
346end
347close(h)
Note: See TracBrowser for help on using the repository browser.

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