Changeset 119 for aedes_readfid.m
 Timestamp:
 Apr 9, 2010, 2:24:32 PM (9 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

aedes_readfid.m
r118 r119 54 54 % 'FastRead' : [{'off'}  'on' ] % Enables an experimental 55 55 % % method for very fast reading 56 % % of fidfiles. This not as56 % % of fidfiles. This is not as 57 57 % % robust as the older 58 58 % % method and can consume a lot 59 59 % % of memory. 60 % 61 % 'FileChunkSize' : [ megabytes ] % Chunk size in megabytes for 62 % % reading fidfiles 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 % 60 71 % 61 72 % 'OutputFile' : filename % Writes the slices into … … 127 138 % % PROCPAR.orient property. 128 139 % 129 % 'RemoveEPIphaseIm' : [ 'on'{'off'}] % Remove the phase image140 % 'RemoveEPIphaseIm' : [{'on'}'off'] % Remove the phase image 130 141 % % from EPI data. This 131 142 % % option only affects EPI 132 143 % % images. 144 % 145 % 'EPIBlockSize' : [integer value] % Block size (number of 146 % % volumes) used for 147 % % processing multireceiver 148 % % EPI data. Default=100 133 149 % 134 150 % Examples: … … 193 209 Dat.FastDataRead = true; 194 210 Dat.precision = 'single'; 195 211 Dat.FileChunkSize = 500; % Chunk size (in MB) for FastRead 196 212 197 213 %% Other defaults … … 204 220 Dat.ReturnRawKspace = false; 205 221 Dat.ReorientEPI = false; 206 Dat.RemoveEPIphaseIm = false; 222 Dat.RemoveEPIphaseIm = true; 223 Dat.EPIBlockSize = 100; 207 224 Dat.OrientImages = true; 208 225 … … 337 354 procpar=varargin{ii+1}; 338 355 339 case 'zeropadding'356 case 'zeropadding' 340 357 if ischar(varargin{ii+1}) 341 358 if strcmpi(varargin{ii+1},'on') … … 352 369 end 353 370 354 case 'phasetable'371 case 'phasetable' 355 372 Dat.phasetable = varargin{ii+1}; 356 373 357 case 'sorting'374 case 'sorting' 358 375 if strcmpi(varargin{ii+1},'on') 359 376 Dat.UsePhaseTable = true; … … 365 382 366 383 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 379 399 380 400 case 'precision' … … 419 439 Dat.RemoveEPIphaseIm = true; 420 440 end 441 case 'epiblocksize' 442 Dat.EPIBlockSize = round(varargin{ii+1}); 421 443 case 'orientimages' 422 444 if strcmpi(varargin{ii+1},'off') … … 1071 1093 d=dir(fopen(file_fid)); 1072 1094 file_sz = d.bytes/1024/1024; % File size in MB 1073 if file_sz< 5001095 if file_sz<Dat.FileChunkSize 1074 1096 nBlocks = 1; 1075 1097 else 1076 nBlocks = ceil(file_sz/ 500); % Read data in 500 MB blocks1098 nBlocks = ceil(file_sz/Dat.FileChunkSize); % Read data in 500 MB blocks 1077 1099 end 1078 1100 1079 1101 % Initialize waitbar 1080 1102 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 1083 1111 end 1084 1112 … … 1132 1160 szh = nVals+hdr.FileHeader.np*hdr.FileHeader.ntraces; 1133 1161 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 1134 1167 1135 1168 % Read the whole data including block headers etc... … … 1161 1194 kspace = reshape(kspace,hdr.FileHeader.np/2,... 1162 1195 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 data1170 % kspace(1:nVals,:)=[];1171 % %kspace=kspace(nVals+1:end,:);1172 %1173 % % Transform to complex values1174 % kspace = reshape(kspace,hdr.FileHeader.np,...1175 % hdr.FileHeader.ntraces,hdr.FileHeader.nblocks);1176 %1177 % % Do DCcorrection if necessary1178 % if ~Dat.DCcorrection  ( nt(1)>1 )1179 % kspace=complex(kspace(1:2:end,:,:),kspace(2:2:end,:,:));1180 % else1181 % kspace=complex(kspace(1:2:end,:,:)hdr.BlockHeaders.lvl,...1182 % kspace(2:2:end,:,:)hdr.BlockHeaders.tlt);1183 % end1184 1196 1185 1197 %% Store and order kspace values … … 1197 1209 1198 1210 if Dat.ShowWaitbar 1199 delete(wb_h) 1211 delete(wb_h) 1212 pause(0.1) 1200 1213 end 1201 1214 … … 1335 1348 % and calculate sumofsquares image 1336 1349 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 1337 1355 nVols = size(kspace,3)/nRcv1; 1338 1356 data = zeros(procpar.nv,procpar.np/2,procpar.ns,nVols+1,'single'); 1339 1357 kssz=size(kspace); 1340 blksz = 500; % Process EPI data in 500 volume blocks1341 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)/nRcvnRef)/blksz); 1342 1360 lnum = length(num2str(nBlocks)); 1343 1361 lnumstr = num2str(lnum); … … 1350 1368 tmp_data = []; 1351 1369 for kk=1:nRcv 1352 inds = [kk ((ii1)*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 ((ii1)*blksz*nRcv+nRcv*nRef+kk):nRcv:min((nRcv*ii*blksz+kk),kssz(3))]; 1357 1371 tmp_kspace = l_ReconstructKspace(kspace(:,:,inds),procpar,Dat); 1358 1372 tmp_data(:,:,:,:,kk) = fftshift(fftshift(fft(fft(tmp_kspace,[],1),[],2),1),2); … … 1361 1375 data(:,:,:,1:size(tmp_data,4)) = sqrt(sum(tmp_data.*conj(tmp_data),5)); 1362 1376 elseif ii==nBlocks 1363 data(:,:,:,((ii1)*blksz+1):(nVols+1)) = sqrt(sum(tmp_data(:,:,:,2:end,:).*conj(tmp_data(:,:,:,2:end,:)),5)); 1377 tmp_data(:,:,:,1:nRef,:)=[]; 1378 data(:,:,:,((ii1)*blksz+1+nRef):(nVols+1)) = sqrt(sum(tmp_data.*conj(tmp_data),5)); 1364 1379 else 1365 data(:,:,:,((ii1)*blksz+1):(ii*blksz)) = sqrt(sum(tmp_data(:,:,:,2:end,:).*conj(tmp_data(:,:,:,2:end,:)),5)); 1380 tmp_data(:,:,:,1:nRef,:)=[]; 1381 data(:,:,:,((ii1)*blksz+1+nRef):(ii*blksz+nRef)) = sqrt(sum(tmp_data.*conj(tmp_data),5)); 1366 1382 end 1367 1383 end … … 1392 1408 % Remove reference image if requested 1393 1409 if Dat.isEPIdata && Dat.RemoveEPIphaseIm 1394 data = data(:,:,:,2:end); 1410 %data = data(:,:,:,2:end); 1411 data(:,:,:,1:nRef)=[]; 1395 1412 end 1396 1413 … … 1413 1430 for ii=1:nRcv 1414 1431 tmp_kspace = l_ReconstructKspace(kspace(:,:,ii),procpar,Dat); 1415 kspace2(:,:,:, ii)=tmp_kspace;1432 kspace2(:,:,:,:,ii)=tmp_kspace; 1416 1433 if Dat.ReturnFTData 1417 1434 tmp_data = l_CalculateFFT(tmp_kspace,procpar,Dat); 1418 data(:,:,:, ii) = tmp_data;1435 data(:,:,:,:,ii) = tmp_data; 1419 1436 end 1420 1437 end … … 1557 1574 %kspace = flipdim(kspace,1); 1558 1575 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 1566 1603 end 1567 1604 end
Note: See TracChangeset
for help on using the changeset viewer.