Multi-frequency reconstruction


Jump to: navigation, search

(Please feel free to edit this page with any notes / tips that you may have, in order to help out your fellow users!)


[edit] Getting the code

The code for multi-frequency reconstruction can be obtained as follows:

 svn checkout

In this code repository, relevant functions are:

  • greconsmultifreq - a grecons wrapper that smoothes fieldmaps and performs multi-frequency reconstruction
  • greconsevaluateb0 - a grecons wrapper that performs many B0 reconstructions

Note that you will also need the grecons executable (contact Gary Glover).

[edit] Information

This is a presentation that may be informative: Kendrick's Vision Lunch Presentation (fixed version)

Here is the beginning of a white paper.

Some advice on smoothing:

[edit] Temporal smoothing

In general, we advise experimenters to ask their subjects to hold their breath for the first few frames of each scan. The reason is that the field maps are derived by comparing frame 1 to frame 2, and if the subject at a different point in the respiratory cycle in these two frames, then the field map will likely be inaccurate.

If your subject held his breath during the beginning of each scan, then we have more confidence in each field map estimate, and we don't need to do as much smoothing of field maps across scans. If the subject did not hold his breath, then each field map measure is likely to be inaccurate, and we need to do more averaging over scans.

If you think the field maps are good, then you probably want to use a short temporal smoothing value (like 2 scans). If you think the field maps are bad, then you might need a longer temporal smoothing value (e.g., 6 or 7 or just 'inf').

[edit] Spatial smoothing

You can determine the best spatial smoothing scale by running several multifrequency recons, each with different smoothing values, and comparing the results. In our experience, spatial smoothing of 7.5 mm works well, and we now just use this value without sweeping out a range of parameters.

[edit] Tutorial

Here is the sample dataset.

 % download the sampledata.tar file and untar it.

 % basic information about this dataset:
 % - 64x64 matrix size, 19 slices, 2.5 mm voxel size (isotropic), spiral-out, 380 volumes.
 % - in-plane anatomicals were acquired with 20 slices (first 19 slices were the ones used for the functional data)
 matrixSz = 64;
 nslices  = 19;
 x        = 2.5; % size of voxel, within plane (mm)
 y        = 2.5; % size of voxel, within plane (mm)
 z        = 2.5; % slice thickness (mm)
 % inspect anatomy
 anat = double(multifile('*.dcm',@dicomread,3));  % load in anatomy
 anat = processmulti(@imresize,anat,[matrixSz matrixSz]);     % downsample
 imwrite(uint8(255*makeimagestack(flipdim(rotatematrix(anat,1,2,1),1),1)),sprintf('anatomy.png'));  % rotate, flip, write to file
 % prep
 grecons = 'grecons14_rev95';  % the grecons utility to call
                               % Note that grecons14 is written specifically for Lucas Center 3T-1 scanner. 
                               % Files from 3T-2 need a different utility (beginning with grecons20)
 files = matchfiles('*.7');    % what files to process
 %  A set of default params
 % If you want a function call with no interactive or iterative steps, you might try this. 
 % Just run it, and you are done. We make no guarantees about the data quality of course. 
 % But these parameters often work for us. If you want to explore parameters and optimize 
 % your recon, then continue to the steps below (and skip this one).
 smoothinsSzInTime = 2;     % (2 scans) We assume the maps of consecutive scans are similar, but far apart scans may not be.
 smoothinsSzInmm   = 7.5;   % (7.5 mm) From experience, this seems to work well. 
 llbands           = makellbands(smoothinsSzInmm, x,y,z, smoothinsSzInTime);
 scans     =[];             % set scans to empty matrix to indicate we want to do all scans
 whichvols =[];             % set whichvols to empty matrix to indicate we want to do all frames
 greconsmultifreq(grecons,files,'final',matrixSz,nslices,llbands,0,[Inf 5],whichvols,scans,1/2);
 % reconstruct functional data using original method 
 %  This step is usually unnecessary. It should produce the same "P.mag" 
 %  files that are produced by default at the scanner. So why would you ever do this? 
 %  1. Maybe you didn't recon at the scanner and don't have .mag files
 %  2. Sanity check to make sure code works 
 %  3. Get some images from the default recon in the same format as the images 
 %     from the multifrequency recon, so that you can easily compare them (and 
 %     feel happy if you improved your images, or sad if you did not)
 % reconstruct functional data using constant reconstructions (ignore fieldmap information)
 %  This step is usually unnecessary. It would be useful in the rare 
 %  case that your field map estimates are corrupted or otherwise
 %  untrustable. In that rare case, the best you can do is to recon
 %  at many different uniform B0 values and choose the B0 value that 
 %  produces the best images. Using this method means giving up on  
 %  correcting for spatial distortions in the B0 field.
 dcrange = -100:10:100;% Choose a range of dc values for constant reconstructions.
                       % One way to estimate an appropriate range is to look
                       % at the parameter1 estimates spit out from the
                       % original reconstrcution. To be safe, choose a range
                       % that is bigger than the range of values in this
                       % plot. 
 whichvols = 6;        % Because we are going to do a separate reconstruction 
                       % for each DC value, we probably do not want to
                       % reconstruct all volumes (this will take a long
                       % time.) We also want to go beyond the first slice
                       % to ensure that the longitudingal relaxation has
                       % reached equilibrium. So we choose a pretty low value
                       % greater than 1. (but why 6? just arbitrary.) 
 % reconstruct functional data using new multi-frequency method:
 %   evaluate different bandwidths (only first 6 volumes)
 %  The purpose of this step is to choose a smoothing kernel for the field
 %  map estimates. To do so we do several multifrequency reconstructions. 
 %  Each recon uses a field map smoothed with a different size smoothing kernel.  
 %  We don't want to waste time reconning every frame of every scan - we will do
 %  that in the next step after we have chosen our parameters. So here we will 
 %  recon only the first few frames and only the first scan, and then look at the 
 %  results. 
 %  For the sample data set we have only one scan and hence do not worry
 %  about smoothing over time. If we have multiple scans, we have to think about
 %  how to smooth over time. This is complicated. If you are confident that there
 %  was little motion across runs and no additional shims, then "inf" is probably
 %  best for temporal smoothing. If there was additional shimming, then scans 
 %  should be reconned completely separately, because the field maps will be different.
 %  If there was significant motion between scans, then you've got some work to do to
 %  figure out the best way to smooth over time, and you are on your own. 
 %  When we have only one scan, 'tweaks' may be important. Tweaks allow us to use the
 %  field map to find the spatial pattern of distortion, and then add an arbitrary DC 
 %  value to the estimated field map. We will look at each tweak value after reconning
 %  and decide which is best. If the field map estimate is good (requires no tweaking),
 %  then the middle tweak value (0, which is the 6th value in our range here) would be best. 
 %  If we have many scans, we probably do not need to tweak since the smoothing over
 %  time should get us the right DC value. 
 smoothingSzInMM = 15 * 2.^([-2:3 Inf]);                 % choose 6 different smoothing bandwiths (in mm)  
 llbands         = makellbands(smoothingSzInMM, x,y,z)  % different bandwidths to try
 tweakrange      = -5:5;                                 % range, in Hz, of DC signal to add to fitted field map
 scans           = 1;                                    % recon only scan 1 (we only have one scan here, but
                                                         % even if we had more, we would do only one to save time 
 for p=size(llbands,1):-1:1
   greconsmultifreq(grecons,files,sprintf('new%02d',p),matrixSz,nslices,llbands(p,:),0,[Inf 5],whichvols,scans,1/2,[],[],[],[],tweakrange);
 % reconstruct functional data using new multi-frequency method:
 %   pick a bandwidth and a tweak value, and then do it for real (all volumes)
 p         = 3;             % this assumes we looked at the images and decided the 3rd level of smoothing had the nicest images
 tw        = tweakrange(4); % this assumes we looked at the images and decided the 4th level of tweaks had the nicest images
 scans     =[];             % set scans to empty matrix to indicate we want to do all scans
 whichvols =[];             % set whichvols to empty matrix to indicate we want to do all frames
 greconsmultifreq(grecons,files,'final',matrixSz,nslices,llbands(p,:),0,[Inf 5],whichvols,scans,1/2,[],[],[],[],tw);
 % view data as a movie
 vol = [];
 for p=1:length(files)
   vol = cat(4,vol,double(permute(loadbinary([files{p} '.mag'],'int16',[matrixSz matrixSz 0 nslices]),[1 2 4 3])));
Personal tools