Changeset 177


Ignore:
Timestamp:
Sep 15, 2011, 10:13:39 AM (8 years ago)
Author:
tjniskan
Message:
  • Added support for reading complex data from bruker 2dseq files.
  • Fixed scaling of 3D data from 2dseq files.
  • Aedes now prompts for handling of complex data (absolute, real

or imag parts) if only complex data is read from a file.

  • Added an option for EPI phase correction in aedes_readvnmr and epi_recon

M aedes_readbruker.m
M aedes_readvnmr.m
M aedes.m
M aedes_revision.m
M vnmr_recon/epi_recon.m

Files:
5 edited

Legend:

Unmodified
Added
Removed
  • aedes.m

    r171 r177  
    175175        end
    176176               
    177         DATA{1}.FTDATA = abs(DATA{1}.KSPACE);
     177        %DATA{1}.FTDATA = abs(DATA{1}.KSPACE);
    178178      end
    179179    elseif isnumeric(DATA) || islogical(DATA)
     
    22922292
    22932293% Check if data file contain only kspace information
    2294 if isfield(DATA{1},'KSPACE') && ...
    2295         (isempty(DATA{1}.FTDATA) & ~isempty(DATA{1}.KSPACE) )
    2296   DATA{1}.FTDATA = abs(DATA{1}.KSPACE);
     2294if length(DATA)==1 && isfield(DATA{1},'KSPACE') && ...
     2295                (isempty(DATA{1}.FTDATA) & ~isempty(DATA{1}.KSPACE) )
     2296        drawnow;
     2297        resp=questdlg(['The inputted DATA structure contains only k-space',...
     2298                ' information. Do you want to view real, imaginary, or absolute',...
     2299                ' part of the complex k-space?'],...
     2300                'Complex data inputted',...
     2301                'Real','Imaginary','Absolute','Absolute');
     2302        if strcmpi(resp,'Real')
     2303                DATA{1}.FTDATA = real(DATA{1}.KSPACE);
     2304        elseif strcmpi(resp,'Imaginary')
     2305                DATA{1}.FTDATA = imag(DATA{1}.KSPACE);
     2306        elseif strcmpi(resp,'Absolute')
     2307                DATA{1}.FTDATA = abs(DATA{1}.KSPACE);
     2308        else
     2309                DATA=[];
     2310                return
     2311        end
    22972312end
    22982313
  • aedes_readbruker.m

    r176 r177  
    5252%        'wbar'        : [ {'on'} | 'off' ]     % Show/hide waitbar
    5353%
     54%        'Precision'   : [{'double'}|'single']  % The precision for the
     55%                                               % reconstructed and scaled
     56%                                               % data. (default='double')
     57%
    5458%
    5559% Examples:
     
    8084Dat.showWbar = true;
    8185Dat.precision = 'double';
     86Dat.return = 1;
    8287
    8388if nargin == 0 || isempty(filename)
     
    95100if strcmpi(fn,'fid')
    96101        isRawData = true;
     102        data_format = aedes_getdataformat(filename);
     103        if strcmpi(data_format,'vnmr')
     104                error('aedes_readbruker cannot read Varian FID-files. Use aedes_readvnmr instead.');
     105        end
    97106elseif strcmpi(fn,'2dseq')
    98107        isRawData = false;
     
    103112% Check varargin length
    104113if rem(length(varargin),2)==1
    105         error('Arguments may only include filename and param/value pairs.')
     114        error('Input arguments may only include filename and param/value pairs.')
    106115end
    107116
    108117% Parse varargin
    109118for ii=1:2:length(varargin)
    110        
     119        param = varargin{ii};
     120        if ~ischar(param) || isempty(param)
     121                error('Parameters must be strings.')
     122        end
     123        value = varargin{ii+1};
     124        switch lower(param)
     125                case 'wbar'
     126                        Dat.wbar = value;
     127                case 'precision'
     128                        Dat.precision = value;
     129                case 'return'
     130                        if ~isnumeric(value) || ~isscalar(value) || isempty(value) || ...
     131                                        isnan(value) || ~isreal(value) || (value-floor(value)) ~= 0 || ...
     132                                        value < 1 || value > 4
     133                                error('Invalid return value. Return value can be an integer in the range 1..4')
     134                        end
     135                        Dat.return = value;
     136                otherwise
     137                        error('Unknown parameter "%s".',param)
     138        end
    111139end
    112140
     
    127155        pe1_table = [];
    128156end
    129 if isempty(data)
     157if isempty(data) && isempty(kspace)
    130158        DATA = [];
    131159        error(msg);
     
    424452                slope = hdr.reco.RECO_map_slope;
    425453                offset = hdr.reco.RECO_map_offset;
     454                if length(hdr.reco.RECO_size) == 2
     455                        nDim = 2;
     456                elseif length(hdr.reco.RECO_size) == 3
     457                        nDim = 3;
     458                end
    426459                for ii=1:length(slope)
    427                         data(:,:,ii) = data(:,:,ii)/slope(ii)+offset(ii);
     460                        if nDim == 2
     461                                data(:,:,ii) = data(:,:,ii)/slope(ii)+offset(ii);
     462                        else
     463                                data(:,:,:,ii) = data(:,:,:,ii)/slope(ii)+offset(ii);
     464                        end
    428465                end
    429466        end
     
    438475else
    439476        % - Read complex data -----------------------------------
    440         msg = 'Reading complex data from 2dseq has not been implemented yet.';
    441         if Dat.showWbar
    442                 close(wbh);
    443         end
     477        [data,count] = fread(fid,inf,precision);
    444478        fclose(fid);
    445         return
     479       
     480        % Reshape and permute to correct size and orientation
     481        data = reshape(data,[hdr.d3proc.IM_SIX,...
     482                hdr.d3proc.IM_SIY,...
     483                hdr.d3proc.IM_SIZ,...
     484                hdr.d3proc.IM_SIT]);
     485        data = permute(data,[2 1 3 4]);
     486       
     487        % Map integer data to single or double values
     488        if ~strcmpi(hdr.reco.RECO_wordtype,'_32BIT_FLOAT')
     489                slope = hdr.reco.RECO_map_slope;
     490                offset = hdr.reco.RECO_map_offset;
     491                if length(hdr.reco.RECO_size) == 2
     492                        nDim = 2;
     493                elseif length(hdr.reco.RECO_size) == 3
     494                        nDim = 3;
     495                end
     496                for ii=1:length(slope)
     497                        if nDim == 2
     498                                data(:,:,ii) = data(:,:,ii)/slope(ii)+offset(ii);
     499                        else
     500                                data(:,:,:,ii) = data(:,:,:,ii)/slope(ii)+offset(ii);
     501                        end
     502                end
     503        end
     504       
     505        data = reshape(data,hdr.reco.RECO_size(1),...
     506                hdr.reco.RECO_size(2),...
     507                hdr.reco.RECO_size(3),[]);
     508       
     509        kspace = complex(data(:,:,:,1),data(:,:,:,2));
     510        data = [];
    446511end
    447512
  • aedes_readvnmr.m

    r165 r177  
    153153%        'EPIPhasedArrayData' : ['on'|{'off'}]  % Return data from
    154154%                                               % individual receivers from
    155 %                                               % phased array EPI files.
     155%                                               % phased array EPI files.
     156%
     157%        'EPI_PC' : [{'auto'}|'pointwise'|      % Phase correction for EPI.
     158%                     'triple'|'off']           % 'auto' = Choose correction
     159%                                               % based on procpar.
     160%                                               % 'pointwise' = Pointwise
     161%                                               % single reference.
     162%                                               % 'triple' = Use triple
     163%                                               % reference correction
     164%                                               % 'off' = Do not perform
     165%                                               % phase correction.
     166%                                                         
    156167%
    157168% Examples:
     
    226237Dat.EPIBlockSize = 300;
    227238Dat.EPIPhasedArrayData = false;
     239Dat.EPI_PC = 'auto';
    228240Dat.OrientImages = true;
    229241Dat.SumOfSquares = 1;
     
    439451        else
    440452          Dat.EPIPhasedArrayData = false;
    441         end
     453                                end
     454                        case 'epi_pc'
     455                               
     456                                if any(strcmpi(varargin{ii+1},...
     457                                                {'off','auto','triple','pointwise'}))
     458                                        Dat.EPI_PC = varargin{ii+1};
     459                                else
     460                                        if ~Dat.ShowErrors
     461                                                msg_out=sprintf('Unknown EPI phase correction type "%s".',varargin{ii+1});
     462                                        else
     463                                                error('Unknown EPI phase correction type "%s".',varargin{ii+1})
     464                                        end
     465                                        return
     466                                end
    442467      case 'orientimages'
    443468        if strcmpi(varargin{ii+1},'off')
  • aedes_revision.m

    r176 r177  
    9393% bash-script every time it is called so that this file "aedes_revision.m" is
    9494% always in the list of committed files. DO NOT EDIT THE NEXT LINE!!!
    95 % - Svn Hook -
     95% - SVN Hook -
  • vnmr_recon/epi_recon.m

    r170 r177  
    6868        end
    6969       
    70         % Number of reference images
    71         if isfield(procpar,'epiref_type') && ~isempty(procpar.epiref_type) && ...
    72                         strcmpi(procpar.epiref_type,'triple')
     70        % EPI phase correction
     71        doPhaseCorrection = true;
     72        if strcmpi(Dat.EPI_PC,'auto')
     73                % Select phase correction using procpar
     74                if isfield(procpar,'epi_pc')
     75                        switch lower(procpar.epi_pc{1})
     76                                case 'triple_ref'
     77                                        nRef = 3;
     78                                case 'scaled_triple'
     79                                        % Issue a warning because unscaled triple ref is used anyway...
     80                                        warning(['Using unscaled TRIPLE REFERENCE phase correction although "%s" ',...
     81                                                'is indicated in PROCPAR.EPI_PC.'],procpar.epi_pc{1});
     82                                        nRef = 3;
     83                                case {'linear','quadratic',...
     84                                                'center_pair','pairwise','first_pair'}
     85                                        % Issue a warning because pointise correction is used anyway...
     86                                        warning(['Using POINTWISE phase correction although "%s" ',...
     87                                                'is indicated in PROCPAR.EPI_PC.'],procpar.epi_pc{1});
     88                                        nRef = 1;
     89                                case 'pointwise'
     90                                        nRef = 1;
     91                                otherwise
     92                                        warning('Unknown EPI_PC "%s". Skipping phase correction,',procpar.epi_pc{1})
     93                        end
     94                else
     95                        % Don't use phase correction
     96                        warning('Cannot find EPI_PC parameter in PROCPAR. Skipping phase correction...');
     97                        doPhaseCorrection = false;
     98                end
     99        elseif strcmpi(Dat.EPI_PC,'off')
     100                nRef = 0;
     101                doPhaseCorrection = false;
     102        elseif strcmpi(Dat.EPI_PC,'triple')
    73103                nRef = 3;
    74         else
     104        elseif strcmpi(Dat.EPI_PC,'pointwise')
    75105                nRef = 1;
    76106        end
     107       
     108%       % Number of reference images
     109%       if isfield(procpar,'epiref_type') && ~isempty(procpar.epiref_type) && ...
     110%                       strcmpi(procpar.epiref_type,'triple')
     111%               nRef = 3;
     112%       else
     113%               nRef = 1;
     114%       end
    77115       
    78116        % Number of segments
     
    82120                nSeg = 1;
    83121        end
     122               
     123       
    84124        %procpar.image = [1 1 0 -2 -1 1 1 1 1];
    85125        % Get index to first reference image
     
    87127        ref_start = ref_start(1);
    88128        nDummys = 0;
    89         if isempty(ref_start)
     129        if isempty(ref_start) && doPhaseCorrection
    90130                kspace = [];
    91131                data = [];
     
    121161        end
    122162       
     163        % Detect partial fourier EPI
     164        if isfield(procpar,'fract_ky') && procpar.fract_ky~=nv/2
     165                nv_full = nv;
     166                nv = nv/2+procpar.fract_ky;
     167        else
     168                nv_full = nv;
     169        end
     170       
    123171        % Calculate indexes for navigator and data echoes
    124172        if nNav > 0
     
    134182       
    135183        if Dat.EPIPhasedArrayData
    136                 data = zeros(nv,np,ns,nVols,nRcvrs,'single');
    137         else
    138                 data = zeros(nv,np,ns,nVols,'single');
     184                data = zeros(nv_full,np,ns,nVols,nRcvrs,'single');
     185        else
     186                data = zeros(nv_full,np,ns,nVols,'single');
    139187        end
    140188       
     
    148196        kspace = reshape(kspace,np,nv+nNav*nSeg,ns,nRcvrs,nVols+nRef+nDummys);
    149197       
     198       
    150199        ref1_ind = find(procpar.image==0,1);
    151200        ref2_ind = find(procpar.image==-2,1);
    152201        im1_ind = find(procpar.image==-1,1);
    153202       
    154         ref_data = kspace(:,DataInd,:,:,[ref1_ind ref2_ind im1_ind]);
    155         if ~isempty(pe_table)
    156                 ref_data(:,pe_table,:,:,:) = ref_data;
    157         end
    158         % Flip even kspace lines
    159         ref_data(:,2:2:end,:,:,:) = flipdim(ref_data(:,2:2:end,:,:,:),1);
    160         %kspace(:,:,:,:,1:nRef)=[];
    161        
    162         % Calculate phase corrections for all receivers
    163         phase_e = ref_data(:,:,:,:,1);
    164         rev_phase = ref_data(:,:,:,:,1);
    165         for ii=1:nRcvrs
    166                 if nRef==1
    167       %ref_im = ref_data(:,:,:,ii,ref1_ind);
    168                         ref_im = ref_data(:,:,:,ii,1);
    169       phase1 = exp(-sqrt(-1)*angle(fft(ref_im,[],1)));
    170                         rev_phase = 1;
    171                         phase_e(:,:,:,ii,:) = phase1;
    172                 elseif nRef==3
    173                         %ref1_ind = find(procpar.image==0,1);
    174                         %ref1 = ref_data(:,:,:,ii,ref1_ind);
    175                         ref1 = ref_data(:,:,:,ii,1);
    176                         phase1 = exp(-sqrt(-1)*angle(fft(ref1,[],1)));
    177                        
    178                         %ref2_ind = find(procpar.image==-2,1);
    179                         %ref2 = flipdim(ref_data(:,:,:,ii,ref2_ind),1);
    180                         ref2 = flipdim(ref_data(:,:,:,ii,2),1);
    181                         phase2 = exp(-sqrt(-1)*angle(fft(ref2,[],1)));
    182                        
    183                         %im1_ind = find(procpar.image==-1,1);
    184                         %im1 = flipdim(ref_data(:,:,:,ii,im1_ind),1);
    185                         im1 = flipdim(ref_data(:,:,:,ii,3),1);
    186                        
    187                        
    188                         rev_phase(:,:,:,ii,:) = fft(im1,[],1).*phase2;
    189                         phase_e(:,:,:,ii,:) = phase1;
     203        if ( isempty(ref2_ind) | isempty(im1_ind) ) && nRef == 3
     204                warning(['Triple reference phase correction was requested, but ',...
     205                        ' < 3 reference images were found. Skipping phase correction.'])
     206                doPhaseCorrection = false;
     207        end
     208               
     209        if doPhaseCorrection
     210
     211                ref_data = kspace(:,DataInd,:,:,[ref1_ind ref2_ind im1_ind]);
     212                if ~isempty(pe_table)
     213                        ref_data(:,pe_table,:,:,:) = ref_data;
     214                end
     215                % Flip even kspace lines
     216                ref_data(:,2:2:end,:,:,:) = flipdim(ref_data(:,2:2:end,:,:,:),1);
     217                %kspace(:,:,:,:,1:nRef)=[];
     218               
     219                % Calculate phase corrections for all receivers
     220                phase_e = ref_data(:,:,:,:,1);
     221                rev_phase = ref_data(:,:,:,:,1);
     222                for ii=1:nRcvrs
     223                        if nRef==1
     224                                %ref_im = ref_data(:,:,:,ii,ref1_ind);
     225                                ref_im = ref_data(:,:,:,ii,1);
     226                               
     227                                phase1 = exp(-sqrt(-1)*angle(fft(ref_im,[],1)));
     228                                rev_phase = 1;
     229                                phase_e(:,:,:,ii,:) = phase1;
     230                               
     231                                %tmp = fft(ref_im,[],1);
     232                                %tmp_max = repmat(max(tmp),[size(ref_im,1),1,1]);
     233                                %phase1 = conj(tmp./tmp_max);
     234                               
     235                                phase_e(:,:,:,ii,:) = phase1;
     236                        elseif nRef==3
     237                                %ref1_ind = find(procpar.image==0,1);
     238                                %ref1 = ref_data(:,:,:,ii,ref1_ind);
     239                                ref1 = ref_data(:,:,:,ii,1);
     240                                phase1 = exp(-sqrt(-1)*angle(fft(ref1,[],1)));
     241                               
     242                                %ref2_ind = find(procpar.image==-2,1);
     243                                %ref2 = flipdim(ref_data(:,:,:,ii,ref2_ind),1);
     244                                ref2 = flipdim(ref_data(:,:,:,ii,2),1);
     245                                phase2 = exp(-sqrt(-1)*angle(fft(ref2,[],1)));
     246                               
     247                                %im1_ind = find(procpar.image==-1,1);
     248                                %im1 = flipdim(ref_data(:,:,:,ii,im1_ind),1);
     249                                im1 = flipdim(ref_data(:,:,:,ii,3),1);
     250                               
     251                               
     252                                rev_phase(:,:,:,ii,:) = fft(im1,[],1).*phase2;
     253                                phase_e(:,:,:,ii,:) = phase1;
     254                        end
    190255                end
    191256        end
     
    220285                % Flip even lines
    221286                tmp_kspace(:,2:2:end,:,:,:) = flipdim(tmp_kspace(:,2:2:end,:,:,:),1);
    222 
    223                 for kk=1:nRcvrs
    224                         if nRef==3
    225                                 for tt=1:length(vol_inds)
    226                                         tmp_kspace(:,:,:,kk,tt) = ifft(rev_phase(:,:,:,kk,:)+fft(tmp_kspace(:,:,:,kk,tt),[],1).*phase_e(:,:,:,kk,:));
     287               
     288                if doPhaseCorrection
     289                        for kk=1:nRcvrs
     290                                if nRef==3
     291                                        for tt=1:length(vol_inds)
     292                                                tmp_kspace(:,:,:,kk,tt) = ifft(rev_phase(:,:,:,kk,:)+fft(tmp_kspace(:,:,:,kk,tt),[],1).*phase_e(:,:,:,kk,:));
     293                                        end
     294                                else
     295                                        for tt=1:length(vol_inds)
     296                                                tmp_kspace(:,:,:,kk,tt) = ifft(fft(tmp_kspace(:,:,:,kk,tt),[],1).*phase_e(:,:,:,kk,:));
     297                                        end
    227298                                end
    228                         else
    229                                 for tt=1:length(vol_inds)
    230                                         tmp_kspace(:,:,:,kk,tt) = ifft(fft(tmp_kspace(:,:,:,kk,tt),[],1).*phase_e(:,:,:,kk,:));
    231                                 end
    232                         end
    233                         tmp_data = fftshift(fftshift(fft(fft(tmp_kspace,[],1),[],2),1),2);
    234                 end
    235                
     299                        end
     300                end
     301               
     302                if nv~=nv_full
     303                        tmp_kspace(:,nv_full,:,:,:) = single(0);
     304                end
     305                tmp_data = fftshift(fftshift(fft(fft(tmp_kspace,[],1),[],2),1),2);
    236306               
    237307                if Dat.EPIPhasedArrayData
Note: See TracChangeset for help on using the changeset viewer.

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