Changeset 195 for aedes_readbruker.m
 Timestamp:
 Feb 22, 2012, 3:45:49 PM (7 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

aedes_readbruker.m
r194 r195 56 56 % % data. (default='double') 57 57 % 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 % % sumofsquares 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 % 58 78 % NOTE: The support for reconstructing raw bruker data is 59 79 % EXPERIMENTAL and should not normally be used. If you don't need … … 89 109 Dat.return = 1; 90 110 Dat.zerofill = 'auto'; 111 Dat.SumOfSquares = 1; 91 112 92 113 if nargin == 0  isempty(filename) … … 131 152 case 'precision' 132 153 Dat.precision = value; 154 case 'sumofsquares' 155 if isscalar(value)  ~any(value~=[1 2]) 156 error('SumOfSquares property cah be 1 (return sumofsquares) 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 133 171 case 'return' 134 172 if ~isnumeric(value)  ~isscalar(value)  isempty(value)  ... … … 594 632 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 595 633 function [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 596 639 597 640 msg = ''; … … 603 646 return 604 647 end 648 649 % Handle EPI data 650 if 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); 656 end 657 605 658 NI = hdr.acqp.NI; % Number of objects 606 659 NSLICES = hdr.acqp.NSLICES; % Number of slices … … 610 663 order = hdr.acqp.ACQ_obj_order; 611 664 665 % Get number of receivers 666 if isfield(hdr,'method') && isfield(hdr.method,'PVM_EncNReceivers') 667 nRcvrs = hdr.method.PVM_EncNReceivers; 668 else 669 nRcvrs = 1; 670 end 671 612 672 % Get phase table 613 673 usePeTable = true; 614 pe1_table= l_GetPeTable(hdr);615 if isempty(pe1_table) 674 [pe1_table,pe2_table] = l_GetPeTable(hdr); 675 if isempty(pe1_table) && isempty(pe2_table) 616 676 usePeTable = false; 617 677 end 618 678 619 679 % 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); 680 if 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); 686 elseif 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); 692 else 693 kspace = []; 694 msg = sprintf('The number of dimensions "%d" is not supported.',nDims); 695 return 629 696 end 630 697 631 698 % Sort echoes 632 699 if 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 634 706 end 635 707 636 708 % Sort object order 637 kspace(:,:,order+1,:) = kspace; 709 if nDims == 2 710 kspace(:,:,order+1,:,:) = kspace; 711 end 638 712 639 713 % Permute to correct orientation … … 698 772 end 699 773 774 % Get number of receivers 775 nRcvrs = size(kspace,5); 776 700 777 % Do Zero filling if requested 701 778 if any(AcqType == [2 3]) 702 779 switch Dat.zerofill 703 780 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 705 798 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 707 808 case 'off' 708 809 padSize = [size(kspace,1) ... 810 size(kspace,2) ... 811 size(kspace,3)]; 709 812 otherwise 710 813 warning('Unknown zerofill option "%s". Skipping zerofilling...',... 711 814 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 827 end 828 829 % Allocate space for Fourier transformed data 830 if nRcvrs>1 && Dat.SumOfSquares==1 831 data_sz = [padSize,size(kspace,4),size(kspace,5)]; 832 data = zeros(data_sz,Dat.precision); 833 else 834 data = zeros(size(kspace),Dat.precision); 835 end 836 714 837 715 838 if 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)); 719 840 elseif 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)); 721 842 elseif 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)); 844 end 845 846 % Calculate sumofsquares image 847 if nRcvrs>1 && Dat.SumOfSquares==1 848 % Calculate sumofsquares 849 data = sqrt(sum(data(:,:,:,:,1:size(data,5)).^2,5)); 850 data=abs(data); 851 elseif 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 856 end 857 724 858 725 859 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 726 860 % Calculate phase table indexes 727 861 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 728 function pt= l_GetPeTable(hdr)862 function [pt1,pt2] = l_GetPeTable(hdr) 729 863 730 864 if ~isfield(hdr,'acqp')  ~isfield(hdr,'method') … … 732 866 end 733 867 868 pt1 = []; 869 pt2 = []; 870 734 871 if 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); 736 873 elseif 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); 875 end 876 877 if isfield(hdr.method,'PVM_EncSteps2') 878 pt2 = hdr.method.PVM_EncSteps2+(min(hdr.method.PVM_EncSteps2(:))+1); 879 end 880 881 882 883 884 885 886 887 888 889
Note: See TracChangeset
for help on using the changeset viewer.