SurfStatReadVol

Reads volumetric files in MINC, ANALYZE, NIFTI or AFNI format. 
 
 Usage: [ data, vol ] = SurfStatReadVol( filenames [,mask [,step  ...
                        [,dirname [, maxmem ] ] ] ] ). 

 filenames = file name with extension .mnc, .img, .nii or .brik as above 
             (n=1), or n x k cell array of file names. Either the data
             files are 4D, or k>1, but not both. If the data files are 4D
             then k=size of the 4th dimension.  
 mask      = 3D logical array of the same size as the unsampled data,
             1=inside, 0=outside. Default is ones, i.e. all voxels.
 step      = 1 x 3 vector of steps for subsampling the volume; only data
             at indices 1, 1+step, 1+2*step, .... is stored. If a scalar,
           = [step step step], or 3 x 3 matrix of 
           = [startx stepx stopx; starty stepy stopy; startz stepz stopz]. 
 dirname   = name of a directory where you have write permission. This
             is only needed if the data is memory mapped. The data is
             memory mapped only if it exceeds maxmem Mb as 4 byte reals.
             SurfStatReadVol then writes the data as 4 byte reals to a
             temporary file in dirname, then memory maps this file to the
             output data. If dirname does not exist, SurfStatReadVol will
             create it. The default is
           = (filenames{1,1} directory)/SurfStat, so you can ignore this
             parameter if you have write permission in the filenames{1,1} 
             directory.
 maxmem    = memory limit in Mb. The default is 64. 

 If filenames is a cell array of file names then 
 data       = n x v matrix of data if k=1, or n x v x k array of data if 
              k>1, or a memory map of same. Note that the mapped file is
              not deleted after you quit MATLAB.
 vol.lat    = 3D logical array, 1=in, 0=out, clamped to the mask.
 vol.vox    = 1 x 3 vector of voxel sizes in mm of the clamped mask.
 vol.origin = position in mm of the first voxel of the clamped mask.
 vol.parent = header info from the first file.
 If n=1, data is double precision; 
 if n>1, data is single precision.

0001 function [data,vol] = SurfStatReadVol(filenames,mask,step,dirname,maxmem);
0002 
0003 %Reads volumetric files in MINC, ANALYZE, NIFTI or AFNI format.
0004 %
0005 % Usage: [ data, vol ] = SurfStatReadVol( filenames [,mask [,step  ...
0006 %                        [,dirname [, maxmem ] ] ] ] ).
0007 %
0008 % filenames = file name with extension .mnc, .img, .nii or .brik as above
0009 %             (n=1), or n x k cell array of file names. Either the data
0010 %             files are 4D, or k>1, but not both. If the data files are 4D
0011 %             then k=size of the 4th dimension.
0012 % mask      = 3D logical array of the same size as the unsampled data,
0013 %             1=inside, 0=outside. Default is ones, i.e. all voxels.
0014 % step      = 1 x 3 vector of steps for subsampling the volume; only data
0015 %             at indices 1, 1+step, 1+2*step, .... is stored. If a scalar,
0016 %           = [step step step], or 3 x 3 matrix of
0017 %           = [startx stepx stopx; starty stepy stopy; startz stepz stopz].
0018 % dirname   = name of a directory where you have write permission. This
0019 %             is only needed if the data is memory mapped. The data is
0020 %             memory mapped only if it exceeds maxmem Mb as 4 byte reals.
0021 %             SurfStatReadVol then writes the data as 4 byte reals to a
0022 %             temporary file in dirname, then memory maps this file to the
0023 %             output data. If dirname does not exist, SurfStatReadVol will
0024 %             create it. The default is
0025 %           = (filenames{1,1} directory)/SurfStat, so you can ignore this
0026 %             parameter if you have write permission in the filenames{1,1}
0027 %             directory.
0028 % maxmem    = memory limit in Mb. The default is 64.
0029 %
0030 % If filenames is a cell array of file names then
0031 % data       = n x v matrix of data if k=1, or n x v x k array of data if
0032 %              k>1, or a memory map of same. Note that the mapped file is
0033 %              not deleted after you quit MATLAB.
0034 % vol.lat    = 3D logical array, 1=in, 0=out, clamped to the mask.
0035 % vol.vox    = 1 x 3 vector of voxel sizes in mm of the clamped mask.
0036 % vol.origin = position in mm of the first voxel of the clamped mask.
0037 % vol.parent = header info from the first file.
0038 % If n=1, data is double precision;
0039 % if n>1, data is single precision.
0040 
0041 if isstr(filenames)
0042     filenames={filenames};
0043 end
0044 
0045 global MAXMEM
0046 if nargin<5
0047     if isempty(MAXMEM)
0048         maxmem=64;
0049     else
0050         maxmem=MAXMEM;
0051     end
0052 end
0053 maxmem=maxmem*2^20;
0054 
0055 d=SurfStatReadVol1(filenames{1,1},0,0);
0056 n=size(filenames,1);
0057 if length(d.dim)==3 | d.dim(4)==1
0058     d.dim(4)=1;
0059     k=size(filenames,2);
0060 else
0061     k=d.dim(4);
0062     kk=size(filenames,2);
0063     if kk>1 & k>1
0064         error(['Can''t read k=' num2str(kk) '>1 files of 4D data, dims= ' num2str(d.dim)]);
0065     end
0066 end
0067 
0068 if nargin<2 | isempty(mask)
0069     mask=logical(ones(d.dim(1:3)));
0070 end
0071 if nargin<3 | isempty(step)
0072     step=1;
0073 end
0074 if length(step)==1
0075     step=[step step step];
0076 end
0077 if isnumeric(step)
0078     i0=1:step(1):d.dim(1);
0079     j0=1:step(2):d.dim(2);
0080     k0=1:step(3):d.dim(3);
0081 else
0082     i0=round((step{1}-d.origin(1))/d.vox(1)+1);
0083     j0=round((step{2}-d.origin(2))/d.vox(2)+1);
0084     k0=round((step{3}-d.origin(3))/d.vox(3)+1);
0085     if isempty(i0)
0086         i0=1:d.dim(1);
0087     else
0088         i0=max(min(i0,d.dim(1)),1);
0089     end
0090     if isempty(j0)
0091         j0=1:d.dim(2);
0092     else
0093         j0=max(min(j0,d.dim(2)),1);
0094     end
0095     if isempty(k0)
0096         k0=1:d.dim(3);
0097     else
0098         k0=max(min(k0,d.dim(3)),1);
0099     end
0100 end    
0101 vol.lat=mask(i0,j0,k0);
0102 mi=find(any(any(vol.lat,2),3));
0103 mj=find(any(any(vol.lat,1),3));
0104 mk=find(any(any(vol.lat,1),2));
0105 bbox=[min(mi) max(mi); min(mj) max(mj); min(mk) max(mk)]';
0106 vol.lat=vol.lat(bbox(1,1):bbox(2,1),bbox(1,2):bbox(2,2),bbox(1,3):bbox(2,3));
0107 steps=zeros(1,3);
0108 if length(i0)>1
0109     steps(1)=i0(2)-i0(1);
0110 end
0111 if length(j0)>1
0112     steps(2)=j0(2)-j0(1);
0113 end
0114 if length(k0)>1
0115     steps(3)=k0(2)-k0(1);
0116 end
0117 vol.vox=d.vox(1:3).*steps;
0118 vol.origin=d.origin+(bbox(1,:)-1).*vol.vox;
0119 vol.parent=d;
0120 
0121 slices=k0(1)+(mk-1)*steps(3);
0122 slicedlat=logical(zeros(d.dim(1),d.dim(2),length(slices)));
0123 is=i0(1)+((bbox(1,1):bbox(2,1))-1).*steps(1);
0124 js=j0(1)+((bbox(1,2):bbox(2,2))-1).*steps(2);
0125 ks=mk-bbox(1,3)+1;
0126 slicedlat(is,js,:)=vol.lat(:,:,ks);
0127 
0128 if n>1
0129     fprintf(1,'%s',[num2str(n) ' files to read, % remaining: 100 ']);
0130     n10=floor(n/10);
0131 end
0132 v=sum(vol.lat(:));
0133 isnum=(n*v*k*4<=maxmem);
0134 if isnum
0135     data=zeros(n,v,k,'single');
0136 else
0137     if nargin<4 | isempty(dirname)
0138         [PATHSTR,NAME,EXT]=fileparts(filenames{1});
0139         dirname=fullfile(PATHSTR,'SurfStat');
0140     end
0141     if ~exist(dirname,'dir')
0142         [SUCCESS,MESSAGE,MESSAGEID]=mkdir(dirname);
0143         if ~SUCCESS
0144             error(MESSAGEID,['Tried to make directory ' dirname ' for memory mapping. \n',...
0145                 'Please specify the name of a directory where you have write permission \n',...
0146                 'as the fourth parameter of SurfStatReadVol - see the help :-) Error: \n',...
0147                 MESSAGE]);
0148         end
0149     end
0150     [PATHSTR,NAME,EXT]=fileparts(tempname);
0151     Filename=fullfile(dirname,NAME);
0152     fid=fopen(Filename,'wb');
0153 end
0154 for i=1:n
0155     if n>1 & rem(i,n10)==0
0156         fprintf(1,'%s',[num2str(round(100*(1-i/n))) ' ']);
0157     end
0158     for j=1:k
0159         if d.dim(4)==1
0160             d=SurfStatReadVol1(filenames{i,j},slices,1);
0161         else
0162             d=SurfStatReadVol1(filenames{i},slices,j);
0163         end
0164         data2=d.data(slicedlat)';
0165         data2(isnan(data2))=0;
0166         if isnum
0167             data(i,:,j)=data2;
0168         else
0169             fwrite(fid,data2,'single');
0170         end
0171     end
0172 end
0173 if n>1
0174     fprintf(1,'%s\n','Done');
0175 end
0176 if ~isnum
0177     fclose(fid);
0178     if k==1
0179         data=memmapfile(Filename,'Format',{'single' [v n] 'Data'},'Writable',true)
0180     else
0181         data=memmapfile(Filename,'Format',{'single' [v k n] 'Data'},'Writable',true)
0182     end
0183 end
0184 
0185 return;
0186 end
0187