source: aedes_readctdata.m

Last change on this file was 178, checked in by tjniskan, 6 years ago
  • Added some support for Bruker spectral data.
  • Fixed a small bug in aedes_readctdata.m.

M aedes_readbruker.m
M aedes_readctdata.m
M aedes_revision.m

File size: 13.4 KB
Line 
1function [DATA,msg] = aedes_readctdata(filename,varargin)
2% AEDES_READCTDATA - Read Gamma Medica CT/SPECT data format
3%   
4%
5% Synopsis:
6%        [DATA,msg] = aedes_readctdata(filename,varargin)
7%
8% Description:
9%
10% Examples:
11%
12% See also:
13%        AEDES
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
32DATA=[];
33msg='';
34ShowWbar = true;
35rotation = 0;
36flipping = '';
37
38if nargin==0 || isempty(filename)
39 
40  % Ask for a file
41  [fname,fpath,findex]=uigetfile({'*.hdr;*.HDR','CT/SPECT header files (*.hdr)';...
42                      '*.xxm;*.XXM','Reconstruction parameter file (*.xxm)';...
43                      '*.*','All Files (*.*)'},...
44                                 'Select a file');
45  if isequal(fname,0) || isequal(fpath,0)
46    if nargout==2
47      msg = 'Cancelled';
48    end
49    return
50  end
51else
52  [fpath,fname,fext]=fileparts(filename);
53  fpath = [fpath,filesep];
54  fname = [fname,fext];
55  if strcmpi(fext,'.hdr')
56    findex = 1;
57  elseif strcmpi(fext,'.xxm')
58    findex = 2;
59  end
60end
61
62% Go through varargin
63for ii=1:2:length(varargin)
64  switch lower(varargin{ii})
65   case 'rotate'
66    rotation = varargin{ii+1};
67   case 'flip'
68    flipping = varargin{ii+1};
69  end
70end
71
72if findex==1
73 
74  % Read header file
75  [hdr,msg] = l_ReadHDRFile([fpath,fname]);
76  if ~isempty(msg)
77    DATA=[];
78    return
79  end
80 
81  if ShowWbar
82    % Show aedes_calc_wait
83    [h,txh]=aedes_calc_wait('Reading SPECT/CT data...');
84    drawnow
85  end
86 
87  % Read data
88  [data,msg] = l_ReadImgFile([fpath,hdr.DataFileName],hdr);
89  if ~isempty(msg)
90    DATA=[];
91    return
92  end
93 
94  if iscell(data)
95   
96    DATA={};
97   
98    for ii=1:length(data)
99      DATA{ii}=struct('DataFormat','spect/ct',...
100                      'FTDATA',data{ii},...
101                      'HDR',hdr);
102    end
103  else
104    DATA.FTDATA = data;
105    DATA.HDR = hdr;
106   
107    % Include file information
108    DATA.HDR.fname = fname;
109    DATA.HDR.fpath = fpath;
110  end
111 
112  if ShowWbar
113    delete(h)
114  end
115 
116elseif findex==2
117 
118  % Look for data files in the same directory
119  %D=dir(fpath);
120  %f_names={D(~[D(:).isdir]).name};
121  %slice_ind=strncmpi(fnames,'slice.',6);
122 
123  % Read parameter file
124  [DATA.HDR,msg]=l_ReadXXMFile([fpath,fname]);
125  if ~isempty(msg)
126    DATA=[];
127    return
128  end
129 
130  % Number of slice files
131  nFiles = DATA.HDR.PARTAG_CUBESIZEZ;
132 
133  slice_files = {};
134  for ii=1:nFiles
135    slice_files{ii}=[fpath,'slice.',sprintf('%04i',ii-1)];
136  end
137 
138% $$$   if all(slice_ind==0)
139% $$$     if nargout==2
140% $$$       msg = 'Could not find slice data files!';
141% $$$     else
142% $$$       error('Could not find slice data files!');
143% $$$     end
144% $$$     return
145% $$$   end
146% $$$   slice_files = {f_names{slice_ind}};
147 
148  % Read data files
149  [DATA.FTDATA,msg]=l_ReadSliceData(slice_files,DATA.HDR,ShowWbar);
150  if ~isempty(msg)
151    DATA=[];
152    return
153  end
154 
155  % Include file information
156  DATA.HDR.fname = fname;
157  DATA.HDR.fpath = fpath;
158 
159end
160
161%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
162% Read *.HDR header file
163%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
164function [HDR,msg]=l_ReadHDRFile(filename)
165
166HDR=[];
167msg='';
168
169% Try to open file for reading
170fid = fopen(filename,'r');
171if fid < 0
172  msg = sprintf('Could not open file "%s"',filename);
173  return
174end
175
176% Read parameter file
177C=textscan(fid,'%s','delimiter','\n');
178fclose(fid); % Close file
179C=C{1};
180if ~strncmp(C{1},'!INTERFILE',10)
181  msg = sprintf('File "%s" is not valid SPECT/CT HDR-file',filename);
182  return
183end
184
185try
186  for ii=1:length(C)
187    if strncmpi(C{ii},'!name of data file',18)
188      if strcmpi(C{ii}(end),'=')
189        HDR.DataFileName = '';
190      else
191        ind=findstr(C{ii},'=');
192        HDR.DataFileName = C{ii}(ind+2:end);
193      end
194    elseif strncmpi(C{ii},'!patient ID',11)
195      if strcmpi(C{ii}(end),'=')
196        HDR.PatientID = '';
197      else
198        ind=findstr(C{ii},'=');
199        HDR.PatientID = C{ii}(ind+2:end);
200      end
201    elseif strncmpi(C{ii},'!study ID',9)
202      if strcmpi(C{ii}(end),'=')
203        HDR.StudyID = '';
204      else
205        ind=findstr(C{ii},'=');
206        HDR.StudyID = C{ii}(ind+2:end);
207      end
208    elseif strncmpi(C{ii},'!matrix size',12)
209      if strncmpi(C{ii},'!matrix size [1]',16)
210        ind=findstr(C{ii},'=');
211        HDR.MatrixSize_1 = str2num(C{ii}(ind+2:end));
212      elseif strncmpi(C{ii},'!matrix size [2]',16)
213        ind=findstr(C{ii},'=');
214        HDR.MatrixSize_2 = str2num(C{ii}(ind+2:end));
215      elseif strncmpi(C{ii},'!matrix size [3]',16)
216        ind=findstr(C{ii},'=');
217        HDR.MatrixSize_3 = str2num(C{ii}(ind+2:end));
218      end
219    elseif strncmpi(C{ii},'!number format',14)
220      if strcmpi(C{ii}(end),'=')
221        HDR.NumberFormat = '';
222      else
223        ind=findstr(C{ii},'=');
224        HDR.NumberFormat = C{ii}(ind+2:end);
225      end
226    elseif strncmpi(C{ii},'!number of bytes per pixel',26)
227      ind=findstr(C{ii},'=');
228      HDR.BytesPerPixel = str2num(C{ii}(ind+2:end));
229    elseif strncmpi(C{ii},'!number of images this frame group',34)
230      ind=findstr(C{ii},'=');
231      HDR.ImagesInFrameGroup = str2num(C{ii}(ind+2:end));
232    elseif strncmpi(C{ii},'!number of frame groups',23)
233      ind=findstr(C{ii},'=');
234      HDR.NbrOfFrameGroups = str2num(C{ii}(ind+2:end));
235     
236      %% Parse frame groups
237      done = false;
238      count = ii+1;
239      framecount = 0;
240      while ~done
241        if strncmpi(C{count},'!frame group number',19)
242          framecount = framecount+1;
243          ind=findstr(C{count},'=');
244          HDR.FrameGroups(framecount).FrameGroupNbr = ...
245              str2num(C{count}(ind+2:end));
246        elseif strncmpi(C{count},'!matrix size [1]',16)
247          ind=findstr(C{count},'=');
248          HDR.FrameGroups(framecount).MatrixSize_1 = str2num(C{count}(ind+2:end));
249        elseif strncmpi(C{count},'!matrix size [2]',16)
250          ind=findstr(C{count},'=');
251          HDR.FrameGroups(framecount).MatrixSize_2 = ...
252              str2num(C{count}(ind+2:end));
253        elseif strncmpi(C{count},'!number format',14)
254          if strcmpi(C{count}(end),'=')
255            HDR.FrameGroups(framecount).NumberFormat = '';
256          else
257            ind=findstr(C{count},'=');
258            HDR.FrameGroups(framecount).NumberFormat = C{count}(ind+2:end);
259          end
260        elseif strncmpi(C{count},'!number of bytes per pixel',26)
261          ind=findstr(C{count},'=');
262          HDR.FrameGroups(framecount).BytesPerPixel = ...
263              str2num(C{count}(ind+2:end));
264        elseif strncmpi(C{count},'!number of images this frame group',34)
265          ind=findstr(C{count},'=');
266          HDR.FrameGroups(framecount).ImagesInFrameGroup = ...
267              str2num(C{count}(ind+2:end));
268        elseif strncmpi(C{count},'image duration',14)
269          ind=findstr(C{count},'=');
270          HDR.FrameGroups(framecount).ImageDuration = ...
271              str2num(C{count}(ind+2:end));
272        elseif strncmpi(C{count},'!maximum pixel count in group',29)
273          ind=findstr(C{count},'=');
274          HDR.FrameGroups(framecount).MaxPixelCount = ...
275              str2num(C{count}(ind+2:end));
276         
277          %% Check if we can bail out from the while loop
278          if HDR.NbrOfFrameGroups==framecount
279            done=true;
280          end
281        elseif strncmpi(C{count},'!END OF INTERFILE',17)
282          %% End of file -> break from while loop
283          done=true;
284        end
285        count =count+1;
286      end
287     
288      %% Break away from the FOR-loop
289      break;
290    else
291      continue
292    end
293  end
294catch
295  msg = sprintf('Error occurred while parsing file "%s"',filename)
296end
297
298
299
300%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
301% Read IMG file
302%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
303function [data,msg]=l_ReadImgFile(filename,HDR);
304
305msg='';
306data = [];
307
308% Try to open img file for reading
309fid = fopen(filename,'r');
310if fid < 0
311  msg = sprintf('Could not open file "%s" for reading',filename);
312  return
313end
314
315% Data size for Frame Group Data
316if isfield(HDR,'FrameGroups')
317 
318  %% Check if the slices in frame groups are of different sizes
319  if all([HDR.FrameGroups(:).MatrixSize_1]==...
320         HDR.FrameGroups(1).MatrixSize_1) && ...
321        all([HDR.FrameGroups(:).MatrixSize_2]==...
322            HDR.FrameGroups(1).MatrixSize_2)
323    equalSize = true;
324   
325    % Allocate variable 'data'
326    data=[];
327  else
328    equalSize = false;
329   
330    % Allocate variable 'data'
331    data={};
332  end
333 
334  %% Read frame groups separately
335  for ii=1:HDR.NbrOfFrameGroups
336   
337    xSize = HDR.FrameGroups(ii).MatrixSize_1;
338    ySize = HDR.FrameGroups(ii).MatrixSize_2;
339    zSize = HDR.FrameGroups(ii).ImagesInFrameGroup;
340   
341    % Get data type
342    if strcmpi(HDR.FrameGroups(ii).NumberFormat,'unsigned integer')
343      if HDR.FrameGroups(ii).BytesPerPixel==2
344        DataStr = '*uint16';
345      elseif HDR.FrameGroups(ii).BytesPerPixel==1
346        DataStr = '*uint8';
347      end
348    elseif strcmpi(HDR.FrameGroups(ii).NumberFormat,'signed integer')
349      if HDR.FrameGroups(ii).BytesPerPixel==2
350        DataStr = '*int16';
351      elseif HDR.FrameGroups(ii).BytesPerPixel==1
352        DataStr = '*int8';
353      end
354    end
355   
356    % Read data
357    try
358      data_tmp = fread(fid,xSize*ySize*zSize,DataStr);
359      data_tmp = reshape(data_tmp,[xSize ySize zSize]);
360     
361      % Orient data
362      data_tmp = permute(data_tmp,[2 1 3]);
363      data_tmp = data_tmp(:,size(data_tmp,2):-1:1,:);
364    catch
365      msg = sprintf('Error while reading data from "%s"',filename);
366     
367      % Close file
368      fclose(fid);
369      return
370    end
371   
372    if equalSize
373      if ii==1
374        data(:,:,end:end+(size(data_tmp,3)-1))=data_tmp;
375      else
376        data(:,:,end+1:end+(size(data_tmp,3)))=data_tmp;
377      end
378    else 
379      for kk=1:size(data_tmp,3)
380        data{end+1}=data_tmp(:,:,kk);
381      end
382    end
383  end
384else
385 
386  % Get data size from header
387  xSize = HDR.MatrixSize_1;
388  ySize = HDR.MatrixSize_2;
389  if isfield(HDR,'MatrixSize_3')
390    zSize = HDR.MatrixSize_3;
391  elseif isfield(HDR,'ImagesInFrameGroup')
392    zSize = HDR.ImagesInFrameGroup;
393    if isfield(HDR,'NbrOfFrameGroups')
394      zSize = zSize*HDR.NbrOfFrameGroups;
395    end
396  else
397    zSize = 1;
398  end
399
400  % Get data type
401  if strcmpi(HDR.NumberFormat,'unsigned integer')
402    if HDR.BytesPerPixel==2
403      DataStr = '*uint16';
404    elseif HDR.BytesPerPixel==1
405      DataStr = '*uint8';
406  end
407  elseif strcmpi(HDR.NumberFormat,'signed integer')
408    if HDR.BytesPerPixel==2
409      DataStr = '*int16';
410    elseif HDR.BytesPerPixel==1
411      DataStr = '*int8';
412    end
413  end
414 
415  % Read data
416  try
417    data = fread(fid,xSize*ySize*zSize,DataStr);
418    data = reshape(data,[xSize ySize zSize]);
419   
420    % Orient data
421    data = permute(data,[2 1 3]);
422    data = data(:,size(data,2):-1:1,:);
423  catch
424    msg = sprintf('Error while reading data from "%s"',filename);
425   
426    % Close file
427    fclose(fid);
428    return
429  end
430
431end
432
433% Close file
434fclose(fid);
435
436
437
438%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
439% Read XXM reconstruction parameter file
440%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
441function [HDR,msg]=l_ReadXXMFile(filename)
442
443
444HDR=[];
445msg='';
446
447% Try to open file for reading
448fid = fopen(filename,'r');
449if fid < 0
450  msg = sprintf('Could not open file "%s"',filename);
451  return
452end
453
454% Read parameter file
455C=textscan(fid,'%s','delimiter','\n');
456fclose(fid); % Close file
457C=C{1};
458
459HDR.Tags={};
460
461% Loop over lines in the file
462for ii=1:length(C)
463 
464  if isempty(C{ii})
465    continue
466  end
467 
468  if strncmpi(C{ii},'BPMODETAG',9)
469    HDR.Tags{end+1} = C{ii};
470    continue;
471  end
472 
473  if strncmpi(C{ii},'****',4)
474    if strncmpi(C{ii},'**** R',6)
475      HDR.Date = C{ii};
476    elseif strncmpi(C{ii},'**** T',6)
477      HDR.TotalReconstructionTime = C{ii};
478    end
479    continue
480  end
481 
482  ind=findstr(C{ii},'=');
483  if isempty(ind)
484    continue
485  end
486 
487  % Extract parameter and value
488  param=C{ii}(1:ind-1);
489  value=C{ii}(ind+1:end);
490 
491  % Test if the parameter is a number
492  val_num=str2num(value);
493  if isempty(val_num);
494    val_num=value;
495  end
496 
497  HDR.(param)=val_num;
498end
499
500%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
501% Read Slice Data
502%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
503function [data,msg]=l_ReadSliceData(slice_files,HDR,ShowWbar)
504
505
506msg='';
507% Get data size from header
508xSize = HDR.PARTAG_CUBESIZEX;
509ySize = HDR.PARTAG_CUBESIZEY;
510nFiles = HDR.PARTAG_CUBESIZEZ;
511
512% Allocate space for data
513data = zeros([xSize,ySize,nFiles],'int16');
514
515% Show waitbar
516if ShowWbar
517  hWbar=aedes_wbar(0,['Reading file 1/',num2str(nFiles)]);
518  tmp_h=findall(hWbar,'string',['Reading file 1/',num2str(nFiles)]);
519  set(tmp_h,'interpreter','none')
520end
521
522% Loop over slice files
523for ii=1:length(slice_files)
524 
525  % Update waitbar
526  if ShowWbar
527    aedes_wbar(ii/nFiles,hWbar,...
528                        {['Reading file ',num2str(ii),'/',num2str(nFiles)],...
529                        slice_files{ii}});
530  end
531 
532  fid = fopen(slice_files{ii},'r');
533  if fid < 0
534    data = [];
535    msg = sprintf('Could not open file "%s" for reading!', ...
536                        slice_files{ii});
537                if ShowWbar
538                        close(hWbar)
539                end
540    return
541  end
542 
543  % Read data
544  data(:,:,ii)=fread(fid,[xSize ySize],'*int16');
545  fclose(fid); % Close file
546end
547
548% Permute data to correct orientation
549data = permute(data,[2 1 3]);
550data = data(:,size(data,2):-1:1,:);
551
552if ShowWbar
553  close(hWbar);
554end
Note: See TracBrowser for help on using the repository browser.

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