Changeset 195 for aedes_readbruker.m


Ignore:
Timestamp:
Feb 22, 2012, 3:45:49 PM (7 years ago)
Author:
tjniskan
Message:
  • Added experimental support for multiple receiver data in aedes_readbruker.
  • Added experimental support for 3D data in aedes_readbruker
  • Added experimental support for zeropadding k-space in aedes_readbruker.

M aedes_readbruker.m
M aedes_revision.m

File:
1 edited

Legend:

Unmodified
Added
Removed
  • aedes_readbruker.m

    r194 r195  
    5656%                                               % data. (default='double')
    5757%
     58%        'ZeroPadding' :  ['off'|'on'|{'auto'}] % 'off' = force
     59%                                               % zeropadding off
     60%                                               % 'on' = force
     61%                                               % zeropadding on (force
     62%                                               % images to be square)
     63%                                               % 'auto' = autoselect
     64%                                               % using relative FOV
     65%                                               % dimensions (default) 
     66%
     67%        'SumOfSquares':  [ {1} | 2 ]          % 1=Return only the
     68%                                               % sum-of-squares image
     69%                                               % for multireceiver
     70%                                               % data (default).
     71%                                               % 2=Return individual
     72%                                               % receiver data
     73%                                               % NOTE: This property has
     74%                                               % no effect on single
     75%                                               % receiver data or when
     76%                                               % reading 2DSEQ files.
     77%
    5878%        NOTE: The support for reconstructing raw bruker data is
    5979%        EXPERIMENTAL and should not normally be used. If you don't need
     
    89109Dat.return = 1;
    90110Dat.zerofill = 'auto';
     111Dat.SumOfSquares = 1;
    91112
    92113if nargin == 0 || isempty(filename)
     
    131152                case 'precision'
    132153                        Dat.precision = value;
     154                case 'sumofsquares'
     155                        if isscalar(value) || ~any(value~=[1 2])
     156                                error('SumOfSquares property cah be 1 (return sum-of-squares) or 2 (return individual receiver data).');
     157                        end
     158                        Dat.SumOfSquares = value;
     159                case 'zeropadding'
     160                         if ischar(varargin{ii+1}) && any(strcmpi(value,{'auto','on','off'}))
     161          if strcmpi(varargin{ii+1},'on')
     162            Dat.zerofill = lower(varargin{ii+1}); % on
     163          elseif strcmpi(varargin{ii+1},'off')
     164            Dat.zerofill = lower(varargin{ii+1}); % off
     165                                        else
     166                                                Dat.zerofill = lower(varargin{ii+1}); % auto
     167                                        end
     168                         else
     169                                 error('ZeroPadding property has to be a string ''auto'', ''on'' or ''off''.')
     170                         end
    133171                case 'return'
    134172                        if ~isnumeric(value) || ~isscalar(value) || isempty(value) || ...
     
    594632%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    595633function [kspace,pe1_table,msg] = l_ReconstructKspace(kspace,hdr,Dat)
     634%
     635% NOTE: The reconstruction of raw bruker data is EXPERIMENTAL and may work
     636% correctly.
     637%
     638
    596639
    597640msg = '';
     
    603646        return
    604647end
     648
     649% Handle EPI data
     650if strncmpi(hdr.acqp.PULPROG,'EPI',3)
     651        kspace = [];
     652        msg = 'EPI reconstruction is not supported.';
     653        return
     654        %epi_matrix_size = hdr.method.PVM_Matrix;
     655        %kspace = reshape(kspace,epi_matrix_size(1),epi_matrix_size(2),NSLICES,NR);
     656end
     657
    605658NI = hdr.acqp.NI; % Number of objects
    606659NSLICES = hdr.acqp.NSLICES; % Number of slices
     
    610663order = hdr.acqp.ACQ_obj_order;
    611664
     665% Get number of receivers
     666if isfield(hdr,'method') && isfield(hdr.method,'PVM_EncNReceivers')
     667        nRcvrs = hdr.method.PVM_EncNReceivers;
     668else
     669        nRcvrs = 1;
     670end
     671
    612672% Get phase table
    613673usePeTable = true;
    614 pe1_table = l_GetPeTable(hdr);
    615 if isempty(pe1_table)
     674[pe1_table,pe2_table] = l_GetPeTable(hdr);
     675if isempty(pe1_table) && isempty(pe2_table)
    616676        usePeTable = false;
    617677end
    618678
    619679% Reshape data so that all echoes are in correct planes
    620 kspace = reshape(kspace,im_size(1),...
    621         phase_factor,NI,im_size(2)/phase_factor,NR);
    622 kspace = permute(kspace,[1 2 4 3 5]);
    623 kspace = reshape(kspace,im_size(1),im_size(2),NI,NR);
    624 
    625 % Handle EPI data
    626 if strncmpi(hdr.acqp.PULPROG,'EPI',3)
    627         epi_matrix_size = hdr.method.PVM_Matrix;
    628         kspace = reshape(kspace,epi_matrix_size(1),epi_matrix_size(2),NSLICES,NR);
     680if nDims == 2
     681        % Handle 2D data
     682        kspace = reshape(kspace,im_size(1),...
     683                nRcvrs,phase_factor,NI,im_size(2)/phase_factor,NR);
     684        kspace = permute(kspace,[1 3 5 4 6 2]);
     685        kspace = reshape(kspace,im_size(1),im_size(2),NI,NR,nRcvrs);
     686elseif nDims == 3
     687        % Handle 3D data
     688        kspace = reshape(kspace,im_size(1),...
     689                nRcvrs,phase_factor,im_size(2)/phase_factor,im_size(3),NR);
     690        kspace = permute(kspace,[1 3 4 5 6 2]);
     691        kspace = reshape(kspace,im_size(1),im_size(2),im_size(3),NR,nRcvrs);
     692else
     693        kspace = [];
     694        msg = sprintf('The number of dimensions "%d" is not supported.',nDims);
     695        return
    629696end
    630697
    631698% Sort echoes
    632699if usePeTable
    633         kspace = kspace(:,pe1_table,:,:);
     700        if ~isempty(pe1_table)
     701                kspace = kspace(:,pe1_table,:,:,:);
     702        end
     703        if nDims == 3 && ~isempty(pe2_table)
     704                kspace = kspace(:,:,pe2_table,:,:);
     705        end
    634706end
    635707
    636708% Sort object order
    637 kspace(:,:,order+1,:) = kspace;
     709if nDims == 2
     710        kspace(:,:,order+1,:,:) = kspace;
     711end
    638712
    639713% Permute to correct orientation
     
    698772end
    699773
     774% Get number of receivers
     775nRcvrs = size(kspace,5);
     776
    700777% Do Zero filling if requested
    701778if any(AcqType == [2 3])
    702779        switch Dat.zerofill
    703780                case 'auto'
    704                        
     781                        % Zeropadding is on "auto", i.e. zeropad to FOV
     782                        fov1 = hdr.acqp.ACQ_fov(1);
     783                        fov2 = hdr.acqp.ACQ_fov(2);
     784                        if AcqType == 3
     785                                fov3 = hdr.acqp.ACQ_fov(3);
     786                        end
     787                        if AcqType==2
     788                                % 2D data
     789                                padSize = [hdr.acqp.ACQ_size(1)/2 ...
     790                                        hdr.acqp.ACQ_size(1)/2*(fov2/fov1) ...
     791                                        size(kspace,3)];
     792                        elseif AcqType==3
     793                                % 3D data
     794                                padSize = [hdr.acqp.ACQ_size(1)/2 ...
     795                                        hdr.acqp.ACQ_size(1)/2*(fov2/fov1) ...
     796                                        hdr.acqp.ACQ_size(1)/2*(fov3/fov1)];
     797                        end
    705798                case 'on'
    706                        
     799                        % Zeropad to square
     800                        if AcqType==1
     801                                padSize = hdr.acqp.ACQ_size(1)/2;
     802                        elseif AcqType==2
     803                                padSize = ones(1,2)*hdr.acqp.ACQ_size(1)/2;
     804                                padSize(3) = size(kspace,3);
     805                        else
     806                                padSize = ones(1,3)*hdr.acqp.ACQ_size(1)/2;
     807                        end
    707808                case 'off'
    708                        
     809                        padSize = [size(kspace,1) ...
     810                                size(kspace,2) ...
     811                                size(kspace,3)];
    709812                otherwise
    710813                        warning('Unknown zerofill option "%s". Skipping zerofilling...',...
    711814                                Dat.zerofill);
    712         end
    713 end
     815                        padSize = [size(kspace,1) ...
     816                                size(kspace,2) ...
     817                                size(kspace,3)];
     818        end
     819       
     820        % Pad with zeros
     821        ks_sz = [size(kspace,1) ...
     822                size(kspace,2) ...
     823                size(kspace,3)];
     824        if any(padSize>ks_sz)
     825                kspace(padSize(1),padSize(2),padSize(3))=0;
     826        end
     827end
     828
     829% Allocate space for Fourier transformed data
     830if nRcvrs>1 && Dat.SumOfSquares==1
     831  data_sz = [padSize,size(kspace,4),size(kspace,5)];
     832  data = zeros(data_sz,Dat.precision);
     833else
     834  data = zeros(size(kspace),Dat.precision);
     835end
     836
    714837
    715838if AcqType==1
    716         %data = flipdim(abs(fftshift(fft(kspace,[],1),1)),1);
    717         data = abs(fftshift(fft(flipdim(kspace,1),[],1),1));
    718   %data = circshift(flipdim(abs(fftshift(fft(kspace,[],1),1)),1),[1 0 0]);
     839  data(:,:,:,:,1:size(data,5)) = abs(fftshift(fft(kspace,[],1),1));
    719840elseif AcqType==2
    720   data = abs(fftshift(fftshift(fft(fft(kspace,[],1),[],2),1),2));
     841  data(:,:,:,:,1:size(data,5)) = abs(fftshift(fftshift(fft(fft(kspace,[],1),[],2),1),2));
    721842elseif AcqType==3
    722   data = abs(fftshift(fftshift(fftshift(fft(fft(fft(kspace,[],1),[],2),[],3),1),2),3));
    723 end
     843  data(:,:,:,:,1:size(data,5)) = abs(fftshift(fftshift(fftshift(fft(fft(fft(kspace,[],1),[],2),[],3),1),2),3));
     844end
     845
     846% Calculate sum-of-squares image
     847if nRcvrs>1 && Dat.SumOfSquares==1
     848  % Calculate sum-of-squares
     849  data = sqrt(sum(data(:,:,:,:,1:size(data,5)).^2,5));
     850  data=abs(data);
     851elseif nRcvrs>1 && Dat.SumOfSquares==2
     852        % Move individual receiver data to 4th dimension if possible...
     853        if size(data,4) == 1
     854                data = reshape(data,size(data,1),size(data,2),size(data,3),size(data,5));
     855        end
     856end
     857
    724858
    725859%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    726860% Calculate phase table indexes
    727861%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    728 function pt = l_GetPeTable(hdr)
     862function [pt1,pt2] = l_GetPeTable(hdr)
    729863
    730864if ~isfield(hdr,'acqp') || ~isfield(hdr,'method')
     
    732866end
    733867
     868pt1 = [];
     869pt2 = [];
     870
    734871if isfield(hdr.acqp,'ACQ_spatial_phase_1')
    735         [s,pt] = sort(hdr.acqp.ACQ_spatial_phase_1);
     872        [s,pt1] = sort(hdr.acqp.ACQ_spatial_phase_1);
    736873elseif isfield(hdr.method,'PVM_EncSteps1')
    737         pt = hdr.method.PVM_EncSteps1+(-min(hdr.method.PVM_EncSteps1(:))+1);
    738 else
    739         pt = [];
    740 end
    741 
    742 
    743 
    744 
    745 
    746 
    747 
    748 
    749 
    750 
    751 
     874        pt1 = hdr.method.PVM_EncSteps1+(-min(hdr.method.PVM_EncSteps1(:))+1);
     875end
     876
     877if isfield(hdr.method,'PVM_EncSteps2')
     878        pt2 = hdr.method.PVM_EncSteps2+(-min(hdr.method.PVM_EncSteps2(:))+1);
     879end
     880
     881
     882
     883
     884
     885
     886
     887
     888
     889
Note: See TracChangeset for help on using the changeset viewer.

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