Changeset 119 for aedes_readfid.m


Ignore:
Timestamp:
Apr 9, 2010, 2:24:32 PM (9 years ago)
Author:
tjniskan
Message:
  • Added support for triple-reference EPI data.
  • Added two new properties for aedes_readfid:
    • FileChunkSize? for reducing memory requirements when reading large fid-files (default FileChunkSize?=500 MB)
    • EPIBlockSize for reducing memory requirements when reading large multireceiver EPI (default=100 volumes)
  • Tried to speed up smoothing in fMRI analysis by performing the

operation in frequency space. Still needs some work...

M misclib/fmri_smooth.m
M misclib/fmri_analysis.m
M aedes_readfid.m
M aedes_revision.m

File:
1 edited

Legend:

Unmodified
Added
Removed
  • aedes_readfid.m

    r118 r119  
    5454%        'FastRead'    :  [{'off'} | 'on' ]     % Enables an experimental
    5555%                                               % method for very fast reading
    56 %                                               % of fid-files. This not as
     56%                                               % of fid-files. This is not as
    5757%                                               % robust as the older
    5858%                                               % method and can consume a lot
    5959%                                               % of memory.
     60%
     61%        'FileChunkSize' : [ megabytes ]        % Chunk size in megabytes for
     62%                                               % reading fid-files
     63%                                               % (default=500 MB).
     64%                                               % Lowering the chunk size
     65%                                               % helps to conserve memory
     66%                                               % when reading large files,
     67%                                               % but is slower. This
     68%                                               % property has no effect if
     69%                                               % FastRead='off'.
     70%                                               
    6071%
    6172%        'OutputFile'  :  filename              % Writes the slices into
     
    127138%                                               % PROCPAR.orient property.
    128139%
    129 %        'RemoveEPIphaseIm' : ['on'|{'off'}]    % Remove the phase image
     140%        'RemoveEPIphaseIm' : [{'on'}|'off']    % Remove the phase image
    130141%                                               % from EPI data. This
    131142%                                               % option only affects EPI
    132143%                                               % images.
     144%
     145%        'EPIBlockSize' : [integer value]       % Block size (number of
     146%                                               % volumes) used for
     147%                                               % processing multireceiver
     148%                                               % EPI data. Default=100
    133149%
    134150% Examples:
     
    193209Dat.FastDataRead = true;
    194210Dat.precision = 'single';
    195 
     211Dat.FileChunkSize = 500; % Chunk size (in MB) for FastRead
    196212
    197213%% Other defaults
     
    204220Dat.ReturnRawKspace = false;
    205221Dat.ReorientEPI = false;
    206 Dat.RemoveEPIphaseIm = false;
     222Dat.RemoveEPIphaseIm = true;
     223Dat.EPIBlockSize = 100;
    207224Dat.OrientImages = true;
    208225
     
    337354      procpar=varargin{ii+1};
    338355     
    339      case 'zeropadding'
     356      case 'zeropadding'
    340357      if ischar(varargin{ii+1})
    341358        if strcmpi(varargin{ii+1},'on')
     
    352369      end
    353370     
    354      case 'phasetable'
     371      case 'phasetable'
    355372      Dat.phasetable = varargin{ii+1};
    356373     
    357      case 'sorting'
     374      case 'sorting'
    358375          if strcmpi(varargin{ii+1},'on')
    359376        Dat.UsePhaseTable = true;
     
    365382     
    366383         case 'fastread'
    367            if strcmpi(varargin{ii+1},'on')
    368                  Dat.FastDataRead = true;
    369            else
    370                  Dat.FastDataRead = false;
    371            end
    372            
    373      case 'wbar'
    374       if strcmpi(varargin{ii+1},'on')
    375         Dat.ShowWaitbar = 1;
    376       else
    377         Dat.ShowWaitbar = 0;
    378           end
     384     if strcmpi(varargin{ii+1},'on')
     385       Dat.FastDataRead = true;
     386     else
     387       Dat.FastDataRead = false;
     388     end
     389     
     390      case 'filechunksize'
     391        Dat.FileChunkSize = round(varargin{ii+1});
     392       
     393      case 'wbar'
     394        if strcmpi(varargin{ii+1},'on')
     395          Dat.ShowWaitbar = 1;
     396        else
     397          Dat.ShowWaitbar = 0;
     398        end
    379399         
    380400         case 'precision'
     
    419439          Dat.RemoveEPIphaseIm = true;
    420440        end
     441      case 'epiblocksize'
     442        Dat.EPIBlockSize = round(varargin{ii+1});
    421443      case 'orientimages'
    422444        if strcmpi(varargin{ii+1},'off')
     
    10711093  d=dir(fopen(file_fid));
    10721094  file_sz = d.bytes/1024/1024; % File size in MB
    1073   if file_sz<500
     1095  if file_sz<Dat.FileChunkSize
    10741096    nBlocks = 1;
    10751097  else
    1076     nBlocks = ceil(file_sz/500); % Read data in 500 MB blocks
     1098    nBlocks = ceil(file_sz/Dat.FileChunkSize); % Read data in 500 MB blocks
    10771099  end
    10781100 
    10791101  % Initialize waitbar
    10801102  if Dat.ShowWaitbar
    1081         wb_h = aedes_calc_wait(['Reading ',num2str(Dat.AcqType),...
    1082           'D VNMR data (seqcon: "' procpar.seqcon '")']);
     1103    if nBlocks==1
     1104      wb_h = aedes_calc_wait(['Reading ',num2str(Dat.AcqType),...
     1105        'D VNMR data (seqcon: "' procpar.seqcon '")']);
     1106    else
     1107      wb_h = aedes_wbar(1/nBlocks,{['Reading ',num2str(Dat.AcqType),...
     1108        'D VNMR data (seqcon: "' procpar.seqcon '")'],...
     1109        sprintf('Reading block 0/%d',nBlocks)});
     1110    end
    10831111  end
    10841112 
     
    11321160  szh = nVals+hdr.FileHeader.np*hdr.FileHeader.ntraces;
    11331161  for ii=1:nBlocks
     1162    if nBlocks~=1
     1163      aedes_wbar(ii/nBlocks,wb_h,{['Reading ',num2str(Dat.AcqType),...
     1164        'D VNMR data (seqcon: "' procpar.seqcon '")'],...
     1165        sprintf('Reading block %d/%d',ii,nBlocks)});
     1166    end
    11341167   
    11351168    % Read the whole data including block headers etc...
     
    11611194  kspace = reshape(kspace,hdr.FileHeader.np/2,...
    11621195    hdr.FileHeader.ntraces,hdr.FileHeader.nblocks);
    1163 
    1164  
    1165 %   % Read the whole data including block headers etc...
    1166 %   kspace = fread(file_fid,inf,prec_str);
    1167 %   kspace=reshape(kspace,nVals+hdr.FileHeader.np*hdr.FileHeader.ntraces,[]);
    1168 %   
    1169 %   % Remove block headers from the data
    1170 %   kspace(1:nVals,:)=[];
    1171 %   %kspace=kspace(nVals+1:end,:);
    1172 %   
    1173 %   % Transform to complex values
    1174 %   kspace = reshape(kspace,hdr.FileHeader.np,...
    1175 %       hdr.FileHeader.ntraces,hdr.FileHeader.nblocks);
    1176 %   
    1177 %   % Do DC-correction if necessary
    1178 %   if ~Dat.DCcorrection || ( nt(1)>1 )
    1179 %     kspace=complex(kspace(1:2:end,:,:),kspace(2:2:end,:,:));
    1180 %   else
    1181 %     kspace=complex(kspace(1:2:end,:,:)-hdr.BlockHeaders.lvl,...
    1182 %       kspace(2:2:end,:,:)-hdr.BlockHeaders.tlt);
    1183 %   end
    11841196 
    11851197  %% Store and order k-space values
     
    11971209 
    11981210  if Dat.ShowWaitbar
    1199         delete(wb_h)
     1211    delete(wb_h)
     1212    pause(0.1)
    12001213  end
    12011214 
     
    13351348    % and calculate sum-of-squares image
    13361349    nRcv = length(find(procpar.rcvrs{1}=='y'));
     1350    if ~isfield(procpar,'epiref_type') || strcmpi(procpar.epiref_type,'single')
     1351      nRef = 1;
     1352    elseif strcmpi(procpar.epiref_type,'triple')
     1353      nRef = 3;
     1354    end
    13371355    nVols = size(kspace,3)/nRcv-1;
    13381356    data = zeros(procpar.nv,procpar.np/2,procpar.ns,nVols+1,'single');
    13391357    kssz=size(kspace);
    1340     blksz = 500; % Process EPI data in 500 volume blocks
    1341     nBlocks = ceil((size(kspace,3)/nRcv)/blksz);
     1358    blksz = Dat.EPIBlockSize; % Process EPI data in 100 volume blocks (default)
     1359    nBlocks = ceil((size(kspace,3)/nRcv-nRef)/blksz);
    13421360    lnum = length(num2str(nBlocks));
    13431361    lnumstr = num2str(lnum);
     
    13501368      tmp_data = [];
    13511369      for kk=1:nRcv
    1352         inds = [kk ((ii-1)*blksz*nRcv+kk):nRcv:min((nRcv*ii*blksz),kssz(3))];
    1353         if ii==1
    1354           % For first block, phase image is in inds twice
    1355           inds = inds(2:end);
    1356         end
     1370        inds = [kk:nRcv:nRcv*nRef ((ii-1)*blksz*nRcv+nRcv*nRef+kk):nRcv:min((nRcv*ii*blksz+kk),kssz(3))];
    13571371        tmp_kspace = l_ReconstructKspace(kspace(:,:,inds),procpar,Dat);
    13581372        tmp_data(:,:,:,:,kk) = fftshift(fftshift(fft(fft(tmp_kspace,[],1),[],2),1),2);
     
    13611375        data(:,:,:,1:size(tmp_data,4)) = sqrt(sum(tmp_data.*conj(tmp_data),5));
    13621376      elseif ii==nBlocks
    1363         data(:,:,:,((ii-1)*blksz+1):(nVols+1)) = sqrt(sum(tmp_data(:,:,:,2:end,:).*conj(tmp_data(:,:,:,2:end,:)),5));
     1377        tmp_data(:,:,:,1:nRef,:)=[];
     1378        data(:,:,:,((ii-1)*blksz+1+nRef):(nVols+1)) = sqrt(sum(tmp_data.*conj(tmp_data),5));
    13641379      else
    1365         data(:,:,:,((ii-1)*blksz+1):(ii*blksz)) = sqrt(sum(tmp_data(:,:,:,2:end,:).*conj(tmp_data(:,:,:,2:end,:)),5));
     1380        tmp_data(:,:,:,1:nRef,:)=[];
     1381        data(:,:,:,((ii-1)*blksz+1+nRef):(ii*blksz+nRef)) = sqrt(sum(tmp_data.*conj(tmp_data),5));
    13661382      end
    13671383    end
     
    13921408    % Remove reference image if requested
    13931409    if Dat.isEPIdata && Dat.RemoveEPIphaseIm
    1394       data = data(:,:,:,2:end);
     1410      %data = data(:,:,:,2:end);
     1411      data(:,:,:,1:nRef)=[];
    13951412    end
    13961413   
     
    14131430    for ii=1:nRcv
    14141431      tmp_kspace = l_ReconstructKspace(kspace(:,:,ii),procpar,Dat);
    1415       kspace2(:,:,:,ii)=tmp_kspace;
     1432      kspace2(:,:,:,:,ii)=tmp_kspace;
    14161433      if Dat.ReturnFTData
    14171434        tmp_data = l_CalculateFFT(tmp_kspace,procpar,Dat);
    1418         data(:,:,:,ii) = tmp_data;
     1435        data(:,:,:,:,ii) = tmp_data;
    14191436      end
    14201437    end
     
    15571574    %kspace = flipdim(kspace,1);
    15581575   
    1559     % Get the reference image(s)
    1560     ref_im = kspace(:,:,:,1);
    1561    
    1562     % Do phase correction.
    1563     phase_e = exp(-sqrt(-1)*angle(fft(ref_im,[],1)));
    1564     for kk=2:size(kspace,4)
    1565       kspace(:,:,:,kk) = ifft(fft(kspace(:,:,:,kk),[],1).*phase_e,[],1);
     1576    % EPI phase correction -------------------------
     1577    if ~isfield(procpar,'epiref_type') || strcmpi(procpar.epiref_type,'single')
     1578      % Single reference pointwise phase correction
     1579     
     1580      % Get the reference image
     1581      ref_im = kspace(:,:,:,1);
     1582     
     1583      % Do phase correction.
     1584      phase_e = exp(-sqrt(-1)*angle(fft(ref_im,[],1)));
     1585      for kk=2:size(kspace,4)
     1586        kspace(:,:,:,kk) = ifft(fft(kspace(:,:,:,kk),[],1).*phase_e,[],1);
     1587      end
     1588    elseif strcmpi(procpar.epiref_type,'triple')
     1589      % Triple reference pointwise phase correction
     1590     
     1591      % Get the reference images
     1592      ref1 = kspace(:,:,:,1);
     1593      phase1 = exp(-sqrt(-1)*angle(fft(ref1,[],1)));
     1594      ref2 = flipdim(kspace(:,:,:,3),1);
     1595      phase2 = exp(-sqrt(-1)*angle(fft(ref2,[],1)));
     1596      im1 = flipdim(kspace(:,:,:,2),1);
     1597     
     1598      % Correct phase for reversed read gradient image
     1599      rev_phase = fft(im1,[],1).*phase2;
     1600      for kk=4:size(kspace,4)
     1601        kspace(:,:,:,kk)=ifft(rev_phase+fft(kspace(:,:,:,kk),[],1).*phase1);
     1602      end
    15661603    end
    15671604  end
Note: See TracChangeset for help on using the changeset viewer.

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