Changeset 175 for aedes_readbruker.m


Ignore:
Timestamp:
Sep 2, 2011, 3:10:33 PM (8 years ago)
Author:
tjniskan
Message:
  • Added scaling for Bruker data (reconstructed 2dseq -files).
  • Some progress in reconstructing Bruker raw data (FID-files).

M aedes_readbruker.m
M aedes_revision.m

File:
1 edited

Legend:

Unmodified
Added
Removed
  • aedes_readbruker.m

    r174 r175  
    55%
    66% Synopsis:
    7 %        [DATA,msg]=aedes_readbruker(filename,'PropertyName1',value1,'PropertyName2',value2,...)
    8 %        [DATA,msg]=aedes_readbruker(filename,'header')
     7%        DATA=aedes_readbruker(filename,'PropertyName1',value1,'PropertyName2',value2,...)
     8%        DATA=aedes_readbruker(filename,'header')
    99%
    1010% Description:
     
    1818%          |-> HDR
    1919%                |-> FileHeader   (data file header)
    20 %                |-> BlockHeaders (data block headers, not stored by default)
     20%                |-> BlockHeaders (Varian FID file data block headers)
    2121%                |-> fname        (file name)
    2222%                |-> fpath        (file path)
     
    2727%          |-> PHASETABLE         (phasetable)
    2828%
    29 %        The DATA structure is returned as empty (without HDR, FTDATA,
    30 %        KSPACE, and PROCPAR fields) if an error has occurred during
    31 %        reading.
     29%        The DATA structure is returned as empty (without the above described fields)
     30%        if an error has occurred during reading.
    3231%       
    3332%        The first input argument is either a string containing full path to
    34 %        the FID-file or the header structure HDR. Only the data file header
     33%        the FID/2DSEQ-file or an empty string. Only the data file headers
    3534%        can be read by giving a string 'header' as the second input
    3635%        argument.
     
    3938%        empty. The returned data can be adjusted by using the 'return'
    4039%        property and values 1, 2, or 3 (see below for more information).
     40%        The return option is ignored for 2DSEQ files.
    4141%
    4242%        The supported property-value pairs in AEDES_READBRUKER (property strings
     
    5353%
    5454% Examples:
    55 %        [DATA,msg]=aedes_readvnmr(filename)             % Read image data
    56 %        from 'filename'
     55%        DATA=aedes_readbruker(filename)   % Read image data from 'filename'
    5756%
    5857% See also:
    59 %        AEDES_READBRUKERRAW, AEDES_READBRUKERRECO, AEDES
     58%        AEDES, AEDES_READVNMR
    6059
    6160% This function is a part of Aedes - A graphical tool for analyzing
     
    235234reco_file = [fp,filesep,'reco'];
    236235if exist(reco_file,'file')~=2
    237         hdr = [];
    238         msg = sprintf('Cannot find RECO file in "%s".',fp);
    239         return
    240 end
    241 try
    242         hdr.reco = aedes_readjcamp(reco_file);
    243 catch
    244         hdr = [];
    245         msg = sprintf('Error while reading "%s". The error was "%s".',...
    246                 reco_file,lasterr);
    247         return
    248 end
    249 
    250 
     236        warning('Cannot find RECO file in "%s"',fp);
     237else
     238        try
     239                hdr.reco = aedes_readjcamp(reco_file);
     240        catch
     241                hdr = [];
     242                msg = sprintf('Error while reading "%s". The error was "%s".',...
     243                        reco_file,lasterr);
     244                return
     245        end
     246end
     247
     248% Check for Functional Imaging Tool files
     249if exist([fp,filesep,'fun'],'dir') == 7
     250       
     251        % To be implemented...
     252       
     253end
    251254
    252255return
     
    301304% Read fid file
    302305tmp = fread(fid,inf,precision);
     306fclose(fid);
    303307kspace = complex(tmp(1:2:end),tmp(2:2:end));
    304 kspace = reshape(kspace,[hdr.acqp.ACQ_size(1)/2 hdr.acqp.ACQ_size(2) hdr.acqp.NI]);
     308kspace = reshape(kspace,hdr.acqp.ACQ_size(1)/2,[]);
    305309
    306310% Reconstruct k-space
     
    330334
    331335% Get word type
    332 if strcmpi(hdr.reco.RECO_wordtype,'_8BIT_UNSGN_INT')
    333         precision = ['uint8=>',Dat.precision];
    334 elseif strcmpi(hdr.reco.RECO_wordtype,'_16BIT_SGN_INT')
    335         precision = ['int16=>',Dat.precision];
    336 elseif strcmpi(hdr.reco.RECO_wordtype,'_32BIT_SGN_INT')
    337         precision = ['int32=>',Dat.precision];
    338 elseif strcmpi(hdr.reco.RECO_wordtype,'_32BIT_FLOAT')
    339         precision = ['single=>',Dat.precision];
    340 else
    341         msg = sprintf('Unknown word type "%s".',hdr.reco.RECO_wordtype);
    342         return
     336if isfield(hdr,'reco')
     337        if strcmpi(hdr.reco.RECO_wordtype,'_8BIT_UNSGN_INT')
     338                precision = ['uint8=>',Dat.precision];
     339        elseif strcmpi(hdr.reco.RECO_wordtype,'_16BIT_SGN_INT')
     340                precision = ['int16=>',Dat.precision];
     341        elseif strcmpi(hdr.reco.RECO_wordtype,'_32BIT_SGN_INT')
     342                precision = ['int32=>',Dat.precision];
     343        elseif strcmpi(hdr.reco.RECO_wordtype,'_32BIT_FLOAT')
     344                precision = ['single=>',Dat.precision];
     345        else
     346                msg = sprintf('Unknown word type "%s".',hdr.reco.RECO_wordtype);
     347                return
     348        end
     349else
     350        % Try to guess the word type from file size if reco file was not found...
     351        sz = [hdr.d3proc.IM_SIX,...
     352                hdr.d3proc.IM_SIY,...
     353                hdr.d3proc.IM_SIZ,...
     354                hdr.d3proc.IM_SIT];
     355        d=dir(filename);
     356        nBytes = d.bytes;
     357        bytesPerPixel = nBytes/prod(sz);
     358        if bytesPerPixel == 1
     359                precision = ['uint8=>',Dat.precision];
     360        elseif bytesPerPixel == 2
     361                precision = ['int16=>',Dat.precision];
     362        else
     363                precision = ['int32=>',Dat.precision];
     364        end
    343365end
    344366
    345367% Get byteorder
    346 if strcmpi(hdr.reco.RECO_byte_order,'littleEndian')
     368if isfield(hdr,'reco')
     369        if strcmpi(hdr.reco.RECO_byte_order,'littleEndian')
     370                byteorder = 'ieee-le';
     371        else
     372                byteorder = 'ieee-be';
     373        end
     374else
     375        % Use little endian if no information is available
    347376        byteorder = 'ieee-le';
    348 else
    349         byteorder = 'ieee-be';
    350377end
    351378
    352379% Get image type
    353 if strcmpi(hdr.reco.RECO_image_type,'COMPLEX_IMAGE')
    354         isDataComplex = true;
    355 elseif strcmpi(hdr.reco.RECO_image_type,'MAGNITUDE_IMAGE')
     380if isfield(hdr,'reco')
     381        if strcmpi(hdr.reco.RECO_image_type,'COMPLEX_IMAGE')
     382                isDataComplex = true;
     383        elseif strcmpi(hdr.reco.RECO_image_type,'MAGNITUDE_IMAGE')
     384                isDataComplex = false;
     385        else
     386                msg = sprintf('Unknown image type "%s".',hdr.reco.RECO_image_type);
     387                return
     388        end
     389else
     390        % Assume magnitude image if information is not available
    356391        isDataComplex = false;
    357 else
    358         msg = sprintf('Unknown image type "%s".',hdr.reco.RECO_image_type);
    359         return
    360392end
    361393
     
    373405        fclose(fid);
    374406       
    375         % Reshape to correct size
     407        % Reshape and permute to correct size and orientation
    376408        data = reshape(data,[hdr.d3proc.IM_SIX,...
    377409                hdr.d3proc.IM_SIY,...
     
    380412        data = permute(data,[2 1 3 4]);
    381413       
     414        % Map integer data to single or double values
     415        if ~strcmpi(hdr.reco.RECO_wordtype,'_32BIT_FLOAT')
     416                slope = hdr.reco.RECO_map_slope;
     417                offset = hdr.reco.RECO_map_offset;
     418                for ii=1:length(slope)
     419                        data(:,:,ii) = data(:,:,ii)/slope(ii)+offset(ii);
     420                end
     421        end
     422       
    382423else
    383424        % - Read complex data -----------------------------------
    384425        msg = 'Reading complex data from 2dseq has not been implemented yet.';
     426        fclose(fid);
    385427        return
    386428end
     
    395437msg = '';
    396438
     439
     440NI = hdr.acqp.NI; % Number of objects
     441NSLICES = hdr.acqp.NSLICES; % Number of slices
     442NR = hdr.acqp.NR; % Number of repetitions
     443phase_factor = hdr.acqp.ACQ_phase_factor; % scans belonging to a single image
     444nDims = hdr.acqp.ACQ_dim;
     445im_size = hdr.acqp.ACQ_size;
     446order = hdr.acqp.ACQ_obj_order;
     447
    397448% Get phase table
    398 pe_table = hdr.method.PVM_EncSteps1+(-min(hdr.method.PVM_EncSteps1(:))+1);
     449%pe_table = hdr.method.PVM_EncSteps1+(-min(hdr.method.PVM_EncSteps1(:))+1);
     450[s,pe1_table] = sort(hdr.acqp.ACQ_spatial_phase_1);
     451
     452% Reshape data so that all echoes are in correct planes
     453kspace = reshape(kspace,im_size(1)/2,...
     454        phase_factor,NI,length(pe1_table)/phase_factor,NR);
     455kspace = permute(kspace,[1 2 4 3 5]);
     456kspace = reshape(kspace,im_size(1)/2,im_size(2),NI,NR);
    399457
    400458% Sort echoes
    401 kspace(:,pe_table,:,:) = kspace;
     459kspace = kspace(:,pe1_table,:,:);
     460
     461% Sort object order
     462kspace(:,:,order+1,:) = kspace;
     463
     464% Permute to correct orientation
     465kspace = flipdim(flipdim(kspace,1),2);
     466
     467%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     468% Check that the required reconstruction parameters are found in acqp
     469%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     470function [ok,msg] = l_CheckAcqpParams(hdr)
     471
     472ok = false;
     473msg = '';
     474
     475% Minimal list of ACQP parameters for data reconstruction
     476minParamList = {...
     477        'ACQ_dim',...
     478        'ACQ_size',...
     479        'NI',...
     480        'NR',...
     481        'BYTORDA',...
     482        'ACQ_word_size',...
     483        'ACQ_scan_size',...
     484        'ACQ_scan_shift',...
     485        'ACQ_phase_factor',...
     486        'ACQ_rare_factor',...
     487        'ACQ_phase_encoding_mode',...
     488        'ACQ_phase_enc_start',...
     489        'ACQ_spatial_size_1',...
     490        'ACQ_spatial_size_2',...
     491        'ACQ_spatial_phase_1',...
     492        'ACQ_spatial_phase_2',...
     493        'ACQ_obj_order',...
     494        'ACQ_mod',...
     495        'DSPFVS',...
     496        'DSPFIRM',...
     497        'DECIM'};
     498
     499param = fieldnames(hdr.acqp);
     500
     501
    402502
    403503%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
    424524
    425525if AcqType==1
    426   data = abs(fftshift(ifft(kspace,[],1),1));
     526  data = abs(fftshift(fft(kspace,[],1),1));
    427527elseif AcqType==2
    428   data = abs(fftshift(fftshift(ifft(ifft(kspace,[],1),[],2),1),2));
     528  data = abs(fftshift(fftshift(fft(fft(kspace,[],1),[],2),1),2));
    429529elseif AcqType==3
    430   data = abs(fftshift(fftshift(fftshift(ifft(ifft(ifft(kspace,[],1),[],2),[],3),1),2),3));
    431 end
    432 
    433 
    434 
    435 
    436 
    437 
    438 
    439 
    440 
    441 
    442 
    443 
    444 
    445 
    446 
     530  data = abs(fftshift(fftshift(fftshift(fft(fft(fft(kspace,[],1),[],2),[],3),1),2),3));
     531end
     532
     533
     534
     535
     536
     537
     538
     539
     540
     541
     542
     543
     544
     545
     546
Note: See TracChangeset for help on using the changeset viewer.

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