Changeset 152


Ignore:
Timestamp:
Jan 18, 2011, 1:04:28 PM (8 years ago)
Author:
tjniskan
Message:
  • aedes_readvnmr.m almost works with epi data...
  • The ROI averages of fMRI time series are now shown in one column instead of two in plugins/fmri_plugins/basic_fmri_analysis.m

M aedes_readvnmr.m
M plugins/fmri_plugins/basic_fmri_analysis.m
M aedes_revision.m
M vnmr_recon/epi_recon.m

Files:
4 edited

Legend:

Unmodified
Added
Removed
  • aedes_readvnmr.m

    r148 r152  
    820820
    821821% Allocate space for kspace
    822 kspace = complex(zeros(hdr.FileHeader.np/2,...
    823           hdr.FileHeader.ntraces*hdr.FileHeader.nblocks,...
    824           Dat.precision));
     822% kspace = complex(zeros(hdr.FileHeader.np/2,...
     823%           hdr.FileHeader.ntraces*hdr.FileHeader.nblocks,...
     824%           Dat.precision));
     825
     826kspace = complex(zeros(hdr.FileHeader.np/2*hdr.FileHeader.ntraces,...
     827        hdr.FileHeader.nblocks,Dat.precision));
    825828
    826829
     
    10411044    tmp=reshape(tmp,nVals+hdr.FileHeader.np*hdr.FileHeader.ntraces,[]);
    10421045    tmp(1:nVals,:)=[];
    1043     tmp=reshape(tmp,hdr.FileHeader.np,[]);
     1046    %tmp=reshape(tmp,hdr.FileHeader.np,[]);
    10441047
    10451048    if ii==nBlocks
  • aedes_revision.m

    r151 r152  
    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 -
  • plugins/fmri_plugins/basic_fmri_analysis.m

    r98 r152  
    113113  fprintf(1,'---------------------------------------\n');
    114114  fh=figure;
    115   if length(ROI)<=3
    116     nRows = length(ROI);
    117     nCols=1;
    118   elseif length(ROI)==4
    119     nRows = 2;
    120     nCols=2;
    121   elseif length(ROI)<10
    122     nRows = 3;
    123     nCols = ceil(length(ROI)/nRows);
    124   else
    125     nRows = 4;
    126     nCols = ceil(length(ROI)/nRows);
    127   end
     115        nRows = length(ROI);
     116        nCols = 1;
     117        ax_w = 0.9;
     118        ax_l = (1-ax_w)/2;
     119        ax_gap = 0.01;
     120        ax_h = (1-(length(ROI)+5)*ax_gap)/length(ROI);
     121        bgax = axes('parent',fh,...
     122                'units','normal','position',[0 0 1 1],...
     123                'xlim',[0 1],'ylim',[0 1],'visible','off');
     124%   if length(ROI)<=3
     125%     nRows = length(ROI);
     126%     nCols=1;
     127%   elseif length(ROI)==4
     128%     nRows = 2;
     129%     nCols=2;
     130%   elseif length(ROI)<10
     131%     nRows = 3;
     132%     nCols = ceil(length(ROI)/nRows);
     133%   else
     134%     nRows = 4;
     135%     nCols = ceil(length(ROI)/nRows);
     136%   end
    128137  for kk=1:length(ROI)
    129138   
     
    151160   
    152161    % Plot results
    153     ax=subplot(nRows,nCols,kk,'align','parent',fh);
     162    %ax=subplot(nRows,nCols,kk,'align','parent',fh,'fontsize',8);
     163                ax=axes('parent',fh,'units','normal',...
     164                        'position',[ax_l 1-kk*ax_h-kk*ax_gap ax_w ax_h],...
     165                        'yaxislocation','right','layer','top','box','on');
    154166    line(1:length(z),z_norm,'color','k',...
    155167      'parent',ax);
    156168    line(1:length(z),z_hat_norm,...
    157169      'color','r','linewidth',2,'parent',ax);
    158     title(['Time series for ROI: ',ROI(kk).label])
    159     ylabel(ax,'BOLD-%');
     170                text(ax_l-0.01,1-kk*ax_h-kk*ax_gap+ax_h/2,...
     171                        ['ROI: ',ROI(kk).label],'parent',bgax,...
     172                        'rotation',90,'fontsize',8,...
     173                        'horizontalalign','center',...
     174                        'verticalalign','bottom');
     175    %title(['Time series for ROI: ',ROI(kk).label],'fontsize',8)
     176    ylabel(ax,'BOLD-%','fontsize',8);
    160177    set(ax,'xlim',[0 length(z)],...
    161178      'ylim',[min(z_norm)-min(z_norm)*0.05 ...
    162       max(z_norm)+max(z_norm)*0.05]);
     179      max(z_norm)+max(z_norm)*0.05],'fontsize',8);
    163180    for ii=1:length(new_onset)
    164181      xdata = [new_onset(ii) new_onset(ii)+new_durat(ii) ...
     
    169186        'xdata',xdata,'ydata',ydata,'FaceColor','r',...
    170187        'FaceAlpha',0.3,'LineStyle','none');
    171     end
     188                end
     189                if kk~=length(ROI)
     190                        set(ax,'xticklabel',[])
     191                end
     192                set(ax,'layer','top')
     193                if kk==length(ROI)
     194                        xlabel('Time (scans)','fontsize',8)
     195                end
    172196   
    173197    fprintf(1,'BOLD-%% %s: %.3f\n',ROI(kk).label,bold_prc);
    174   end
     198        end
     199        if length(ROI)>=3
     200                tmp_pos1 = get(0,'screensize');
     201                tmp_pos2 = get(fh,'position');
     202                set(fh,'position',[tmp_pos2(1) 100 tmp_pos2(3)*1.5 tmp_pos1(4)-100]);
     203        end
    175204  fprintf(1,'---------------------------------------\n');
    176205else
  • vnmr_recon/epi_recon.m

    r128 r152  
    66% this code reconstructs
    77if nargin==0
    8   kspace = {'epi_se_rapid_sp3','epi','epi_fMRI'};
     8  kspace = {'epi_se_rapid_sp3','epi','epi_fMRI','epip'};
    99  return
    1010end
     
    1313msg_out = '';
    1414
    15 % EPI data is measured with INOVA system and sorted using CORREPI
     15% Number of receivers
     16if isfield(procpar,'rcvrs') && ~isempty(procpar.rcvrs)
     17        nRcvrs = length(procpar.rcvrs{1}=='y');
     18else
     19        nRcvrs = 1;
     20end
     21
     22
     23
     24% =====================================================
     25% Reconstruct EPI sequences meaured with the old INOVA
     26% =====================================================
    1627if isfield(procpar,'readres') && isfield(procpar,'phaseres')
    1728 
     
    3142  kspace = reshape(kspace,[size(kspace,1) size(kspace,2) ...
    3243    tmp_ns size(kspace,3)/tmp_ns]);
     44       
     45       
     46        % ===================================
     47        % Reconstruct VNMRj and EPIP sequences
     48        % ===================================
    3349else
    34  
     50       
     51        if strncmpi(procpar.seqfil,'epip',4) && isfield(procpar,'nphase') && ...
     52                        isfield(procpar,'nread')
     53                isEPIP = true;
     54        else
     55                isEPIP = false;
     56        end
     57       
     58        % Number of navigator echos
     59        if isfield(procpar,'navigator') && strcmpi(procpar.navigator,'y')
     60                nNav = length(procpar.nav_echo);
     61        else
     62                nNav = 0;
     63        end
     64       
     65        % Number of reference images
     66        if isfield(procpar,'epiref_type') && ~isempty(procpar.epiref_type) && ...
     67                        strcmpi(procpar.epiref_type,'triple')
     68                nRef = 3;
     69        else
     70                nRef = 1;
     71        end
     72       
     73        % Number of volumes
     74        nVols = length(procpar.image)-nRef;
     75       
     76        % Number of slices
     77        ns = procpar.ns;
     78       
     79        % Phase sampling
     80        if isEPIP
     81                nv = procpar.nphase;
     82               
     83                % Read sampling
     84                np = procpar.nread/2;
     85        else
     86                nv = procpar.nv;
     87                np = procpar.np/2;
     88        end
     89       
     90        % Get orientation
     91        orient = procpar.orient{1}
     92       
     93        % Remove navigators
     94        if nNav~=0
     95                nav_data = kspace(1:(np*nNav),:);
     96                kspace(1:(np*nNav),:)=[];
     97        else
     98                nav_data = [];
     99        end
     100       
     101        if Dat.EPIPhasedArrayData
     102                data = zeros(nv,np,ns,nVols,nRcv,'single');
     103        else
     104                data = zeros(nv,np,ns,nVols,'single');
     105        end
     106       
     107        % Reshape kspace and detach reference images
     108        if isEPIP
     109                kspace = reshape(kspace,np,nv,nRcvrs,ns,nVols+nRef);
     110        else
     111                kspace = reshape(kspace,np,nv,nRcvrs,ns,nVols+nRef);
     112        end
     113       
     114        % Flip even kspace lines
     115        if ~Dat.FastDataRead
     116                for tt=2:2:size(kspace,2)
     117                        kspace(:,tt,:,:,:) = flipdim(kspace(:,tt,:,:,:),1);
     118                end
     119        else
     120                kspace(:,2:2:end,:,:,:) = flipdim(kspace(:,2:2:end,:,:,:),1);
     121        end
     122        ref_data = kspace(:,:,:,:,1:nRef);
     123        kspace(:,:,:,:,1:nRef)=[];
     124       
     125        % Calculate phase corrections for all receivers
     126        phase_e = ref_data(:,:,:,:,1);
     127        rev_phase = ref_data(:,:,:,:,1);
     128        for ii=1:nRcvrs
     129                if nRef==1
     130      ref_im = ref_data(:,:,ii,:,1);
     131      phase1 = exp(-sqrt(-1)*angle(fft(ref_im,[],1)));
     132                        rev_phase = 1;
     133                        phase_e(:,:,ii,:,:) = phase1;
     134                elseif nRef==3
     135                        ref1_ind = find(procpar.image==0,1);
     136                        ref1 = ref_data(:,:,ii,:,ref1_ind);
     137                        phase1 = exp(-sqrt(-1)*angle(fft(ref1,[],1)));
     138                       
     139                        ref2_ind = find(procpar.image==-2,1);
     140                        ref2 = flipdim(ref_data(:,:,ii,:,ref2_ind),1);
     141                        phase2 = exp(-sqrt(-1)*angle(fft(ref2,[],1)));
     142                       
     143                        im1_ind = find(procpar.image==-1,1);
     144                        im1 = flipdim(ref_data(:,:,ii,:,im1_ind),1);
     145                       
     146                        rev_phase(:,:,ii,:,:) = fft(im1,[],1).*phase2;
     147                        phase_e(:,:,ii,:,:) = phase1;
     148                end
     149        end
     150       
     151        % Reconstruct in blocks
     152        kssz=size(kspace);
     153        blksz = Dat.EPIBlockSize; % Process EPI data in 100 volume blocks (default)
     154        nBlocks = ceil(nVols/blksz);
     155        lnum = length(num2str(nBlocks));
     156        lnumstr = num2str(lnum);
     157        bsl = lnum*2+1;
     158        fprintf(1,'Processing data in blocks of %d volumes\n',blksz)
     159        fprintf(1,['Processing block...%0',lnumstr,'d/%0',lnumstr,'d'],1,nBlocks);
     160       
     161        for ii=1:nBlocks
     162                fprintf(1,repmat('\b',1,bsl));
     163                fprintf(1,['%0',lnumstr,'d/%0',lnumstr,'d'],ii,nBlocks);
     164                tmp_data = [];
     165                if ii==nBlocks
     166                        vol_inds = (ii-1)*blksz+1:nVols;
     167                else
     168                        vol_inds = (ii-1)*blksz+1:ii*blksz;
     169                end
     170                tmp_kspace = kspace(:,:,:,:,vol_inds);
     171               
     172                for kk=1:nRcvrs
     173                        if nRef==3
     174                                for tt=1:length(vol_inds)
     175                                        tmp_kspace(:,:,kk,:,tt) = ifft(rev_phase(:,:,kk,:,:)+fft(kspace(:,:,kk,:,vol_inds(tt)),[],1).*phase_e(:,:,kk,:,:));
     176                                end
     177                        else
     178                                for tt=1:length(vol_inds)
     179                                        tmp_kspace(:,:,kk,:,tt) = ifft(fft(kspace(:,:,kk,:,vol_inds(tt)),[],1).*phase_e(:,:,kk,:,:));
     180                                end
     181                        end
     182                        tmp_data = fftshift(fftshift(fft(fft(tmp_kspace,[],1),[],2),1),2);
     183                end
     184               
     185               
     186                if Dat.EPIPhasedArrayData
     187                        data_block = abs(tmp_data);
     188                        data_block = permute(data_block,[1 2 4 5 3]);
     189                        % Permute to correct orientation
     190                        if strcmpi(orient,'trans90')
     191                                data_block = permute(data_block,[2 1 3 4 5]);
     192                        end
     193                        data(:,:,:,vol_inds,:) = data_block;
     194                else
     195                        data_block = sqrt(sum(tmp_data.*conj(tmp_data),3));
     196                        % Permute to correct orientation
     197                        if strcmpi(orient,'trans90')
     198                                data_block = permute(data_block,[2 1 3 4 5]);
     199                        end
     200                        data(:,:,:,vol_inds) = squeeze(data_block);
     201                end
     202               
     203        end
     204        fprintf(1,'\n')
     205       
    35206end
     207       
     208
     209
     210
     211
     212
     213
     214
     215
     216
     217
     218
Note: See TracChangeset for help on using the changeset viewer.

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