ncigt-fil
Functions
Parallel Imaging Algorithms

Image reconstruction methods for subsampled data acquired using multiple receiver coils. More...

Functions

int apply_vbc_coefs (DATA_OBJ *input, DATA_OBJ *output, DATA_OBJ *vbc_obj)
 
int compute_root_sum_of_squares (DATA_OBJ *input, DATA_OBJ *output, double *max)
 
int compute_virtual_body_coil (DATA_OBJ *input, DATA_OBJ *output, DATA_OBJ *vbc_obj)
 
int export_grappa_coef (GRAPPA_PARMS **gparm, DATA_OBJ *gexport, int ncoil, int nslc)
 
void generate_b1_filter (int fill_fac, short Ny, short Nz_proc, float *samp_dens, float *filt_b1)
 
int grappa (COMPLEX *I, int *PhaseLines, dataheader *hdr, COMPLEX *O, int *Nlines, int **OutPhaseLines, int YKS, int XKS, float radius, int ndk, int *dks_in)
 
int grappa_calc_parms (dataheader *hdr, COMPLEX *O, int YKS, int *indx0, int XKS, int *indx3, int nacs, int *acsl, int kxnum, COMPLEX **gparm, int fullsqr, int direction)
 
int grappa_calc_recon_coef (COMPLEX *O, int *PhaseLines, dataheader *hdr, GRAPPA_PARMS *gparm, float radius)
 
int grappa_for_cine (short Nx, short N_acq, short Ny, int Ny_proc, short Nz_proc, short Nf, short ncoils, short fill_fac, short fr1, int first_line, float pF, KCOORD *kline_c, float *samp_dens, float *weight_pF, COMPLEX *buf2, COMPLEX *b1, COMPLEX *data, COMPLEX *full)
 
int grappa_recon (dataheader *hdr, COMPLEX *O, int YKS, int *indx0, int XKS, int *indx3, int *bins, COMPLEX *gparm, float radius)
 
int grog_generate_eigd (GRAPPA_PARMS **gparm, int ncoil, int nslc)
 
int grog_grid_2Ddata (DATA_OBJ *smpl_x, DATA_OBJ *smpl_y, DATA_OBJ *kdata, DATA_OBJ *src, DATA_OBJ *new_x, DATA_OBJ *new_y, DATA_OBJ *kdata_out, GRAPPA_PARMS *gparms)
 
int import_grappa_coef (GRAPPA_PARMS ***gparm_out, DATA_OBJ *gimport)
 
int inv_sens_matrix (short ncoils, short npts, float thresh, double **a, double *w, double **v, double *wgt_roe, double *wgt, double **S, double **S_buf, double **S_inv, double **S_inv_buf)
 
int rlsgrappa (COMPLEX *I, int *PhaseLines, dataheader *hdr, COMPLEX *O, int *Nlines, int **OutPhaseLines, int YKS, float radius, int ndk, int *dks_in, float lambda)
 
void self_ref_b1_basic (COMPLEX *kspace_data, float b1snr_thr, short Ny, short Nx, short Nz, short ncoils, short Nf, fil_fft_mngr *mngr, COMPLEX *b1)
 
void self_ref_b1_for_cine (COMPLEX *kspace_data, float sigmak, float b1snr_thr, short N_acq, short Nx, short Ny, short Ny_proc, short Nz_proc, short ncoils, short Nf, short fr1, short fill_fac, short first_line, float *filt_b1, KCOORD *kline_o, KCOORD *kline_e, float *samp_dens, fil_fft_mngr *mngr, COMPLEX *bufsp, COMPLEX *b1)
 
int self_ref_b1_via_geyser (COMPLEX *datain, dataheader *hdr, int *PhaseLines, COMPLEX *ptrW, fil_fft_mngr *mngr, int m)
 
int spacerip_for_cine (short Nx, short N_o, short N_e, short N_c, short Ny, int Ny_proc, short Nz_proc, short Nf, short ncoils, float pF, KCOORD *kline_o, KCOORD *kline_e, KCOORD *kline_c, COMPLEX *b1, COMPLEX *buf, COMPLEX *data, COMPLEX *spacerip, short recon_flag)
 
int sprip3d_conv_pts2y (int nph, int *phz_y, int *phz_z, int M, int N, int *Y)
 
void sprip3d_lsqr_init (int nph, int M, int N, int L)
 
void sprip3d_lsqr_post ()
 
void sprip3d_lsqrsolv (int *Y, int nph, int m, int n, int L, COMPLEX *W, COMPLEX *Sk, COMPLEX **I, float *inlams, int lamsize, int maxits, double *xnormout, double *Rout, int tid)
 
void sprip_cgsolv (int *Y, int nph, int n, int L, COMPLEX *Wx0, COMPLEX *u, COMPLEX *x, COMPLEX lambda, int k_max, int tid)
 
void sprip_cgsr_init (int nph, int N, int L, int threadnum)
 
void sprip_cgsr_post ()
 
void sprip_cgsr_recon (float *phz_encode_list, int P, int N, int M, int L, COMPLEX *W, COMPLEX *Sk, double tau, int maxits, COMPLEX *rho, int tmax)
 
void sprip_lsqr_init (int nph, int N, int L, int threadnum, int lamsize, int maxits)
 
void sprip_lsqr_post ()
 
void sprip_lsqr_recon (float *phz_encode_list, int P, int N, int M, int L, COMPLEX *W, COMPLEX *Sk, float *lams, int numlams, COMPLEX *D, int maxits, COMPLEX **rho, double *xout, double *rout, int tmax, float *usedlams)
 
void sprip_lsqrsolv (int *Y, int nph, int N, int L, COMPLEX *Wx0, COMPLEX *Sk, COMPLEX **I, float *lams, int lamsize, COMPLEX *D, int maxits, double *xnormout, double *Rout, int tid, float *lout)
 
void vdsense (short Nx, short N_acq, short Ny, short Ny_proc, short Nz_proc, short Nf, short ncoils, short fill_fac, short fr1, short first_line, float pF, fil_fft_mngr *mngr, KCOORD *kline_o, KCOORD *kline_e, float *samp_dens, float *weight_pF, COMPLEX *buf, DATATYPE *Uarray, COMPLEX *data, COMPLEX *vdsen_ima)
 
void vdsense_apply (short Nx, short Ny, short Nz, short Nf, short ncoils, COMPLEX *sen_before, DATATYPE *Uarray, COMPLEX *sen_after)
 
void vdsense_prep (short Nx, short Ny, short Nz, short ncoils, float thresh, int accel, XCOORD *alias, COMPLEX *b1, DATATYPE *Uarray)
 

Detailed Description

Image reconstruction methods for subsampled data acquired using multiple receiver coils.

Parallel imaging employs multiple receiver coils, which effectively apply a spatial-domain encoding in addition to the Fourier-domain encoding traditionally employed by MR imaging.

With multiple coils, the MR data acquisition can be modeled by the following signal equation [A]:

\[ s_l(\mathbf{k}) = \int_V \rho(\mathbf{r}) W_l(\mathbf{r}) e^{-j \mathbf{k} \cdot \mathbf{r}} d\mathbf{r} \]

One approch to solve this inverse problem computationally is to first discretize the terms and then concatenate the equations for all coils to form a single linear system of equations.

\[ \mathbf{s} = \mathbf{P} \rho \]

where $ \mathbf{P} $ is described by an estimate of the coil sensitivity maps and the phase-encoding pattern employed. The library provides one approach to estimating coil sensitivity maps, using self-referenced data [D].

Two methods in the NC-IGT fast imaging library employ this approach, SENSE [A] and SPACE RIP[F]. The SENSE functions provide reconstruction of uniformly subsampled data, i.e. where the distance is k-space between each phase-encode line is the same. Variable density SENSE [B] can be employed in those cases where low-frequency and high-frequency k-space are each uniformly sub-sampled, but at different rates. E.g. 2x in the low-frequency range, and 4x in the high-frequency range.

For sampling patterns with a wider variety of k-space distances, that is non-uniform subsampling, one can employ SPACE RIP. Two iterative linear system solvers are provided for SPACE RIP reconstructions: CG [E] and LSQR-Hybrid [I]. At high-acceleration factors, it is often useful to employ regularization.

\[ \left[ \begin{array}{c} \mathbf{s} \\ 0 \end{array} \right] = \left[ \begin{array}{c} \mathbf{P} \\ \lambda \mathbf{I} \end{array} \right] \rho \]

The LSQR-Hybrid method can efficiently find all of the solutions corresponding to a number of regularization parameter values. The "best" image as measured by L-curve methods is then presented as the solution. For a single regularization parameter, CG is slightly faster than LSQR-Hybrid, so this method is provided as well.

An alternate approach to solving the inverse problem is to view the acquisition in k-space as a convolution between two spatial-domain functions. This is the approach taken by GRAPPA [G]. This method does not explicitly employ estimates of the coil sensitivity maps. Rather it maps various output data points together to estimate a convolution kernel combining data from all coils.

Modelling GRAPPA as a matrix operator [L], one can modify the target distance for reconstructed points at will using an Eigen decomposition of the operator [M]. Since Ver 2.0, the library includes functions to compute GRAPPA coefficients suitable for GROG, as well as functions to apply GROG in multi-line radial acquisitions as in PROPELLER [N].

The NC-IGT GRAPPA implementation can reconstruct data acquired using non-uniform sampling patterns, although the visual results typically contain more noise than a corresponding SPACE RIP reconstruction. GRAPPA is the method of choice, however, in cases where aliasing due to a reduced-field-of-view is present in addition to aliasing resulting from subsampling.

It is noted as well that combinations of GRAPPA, VD-SENSE, and SPACE RIP, can often produce better images than any single algorithm alone.

References

  1. Pruessmann KP, Weiger M, Scheidegger MB, Boesiger P. SENSE: Sensitivity encoding for fast MRI. Magn Reson Med 1999; 42(5):952-62.

  2. Madore B. Using UNFOLD to remove artifacts in parallel imaging and in partial-Fourier imaging. Magn Reson Med 2002; 48(3):493-501.

  3. Kellman P, Epstein FH, McVeigh ER. Adaptive sensitivity encoding incorporating temporal filtering (TSENSE). Magn Reson Med 2001; 45(5):846-852.

  4. McKenzie CA, Yeh EN, Ohliger MA, Price MD, Sodickson DK. Self-calibrating parallel imaging with automatic coil sensitivity extraction. Magn Reson Med 2002; 47(3):529-538.

  5. Pruessmann KP, Weiger M, Bornert P, Boesiger P. Advances in sensitivity encoding with arbitrary k-space trajectories. Magn Reson Med 2001; 46(4):638-651.

  6. Kyriakos WE, Panych LP, Kacher DF, Westin CF, Bao SM, Mulkern RV, Jolesz FA. Sensitivity profiles from an array of coils for encoding and reconstruction in parallel (SPACE RIP). Magn Reson Med 2000; 44(2):301-308.

  7. Griswold MA, Jakob PM, Heidemann RM, Nittka M, Jellus V, Wang J, Kiefer B, Haase A. Generalized autocalibrating partially parallel acquisitions (GRAPPA). Magn Reson Med 2002; 47(6):1202-1210.

  8. Hoge WS, Brooks DH, Madore B, Kyriakos WE. A tour of accelerated parallel MR imaging from a linear systems perspective. Concepts in MR 2005; 27A(1):17-37

  9. Hoge WS, Kilmer ME, Haker SJ, Brooks DH, Kyriakos WE. Fast regularized reconstruction of non-uniformly subsampled parallel MRI data. in Proc of 2006 IEEE Intl Symp on Biomedical Imaging (ISBI06). Arlington, VA, USA, 2006; 714-717.

  10. Buehrer M, Boesiger P, and Kozerke S. "Virtual body coil calibration for phased-array imaging." 17th ISMRM Scientific Meeting, page 759, Honolulu, HI, 2009.

  11. Buehrer M, Pruessmann KP, Boesiger P, and Kozerke S. Array compression for MRI with large coil arrays. Magn Reson Med, 57(6):1131-1139, June 2007.

  12. Griswold MA, Blaimer M, Breuer F, Heidemann RM, Mueller M, and Jakob PM. Parallel magnetic resonance imaging using the GRAPPA operator formalism. Magn Reson Med, 54(6):1553-1556, Oct 2005

  13. Seiberlich N, Breuer FA, Blaimer M, Barkauskas K, Jakob PM, and Griswold MA. "Non-Cartesian data reconstruction using GRAPPA operator gridding (GROG)." Magn Reson Med, 58(6):1257-1265, Dec 2007.

  14. Pipe JG. "Motion correction with PROPELLER MRI: Application to head motion and free-breathing cardiac imaging". Magn Reson Med, 42(5):963-969, Nov 1999.

Function Documentation

int apply_vbc_coefs ( DATA_OBJ *  input,
DATA_OBJ *  output,
DATA_OBJ *  vbc_obj 
)

apply a set of coefficients calculated by compute_virtual_body_coil to a new set of input data.

If the output object is that same as the input, the input data will be overwritten, and the size adjusted accordingly

Parameters
inputthe input coil data, sized [x-y-c] or [x-y-z-c]
vbc_objreturns the coil combination coeficients, if non-NULL
int compute_root_sum_of_squares ( DATA_OBJ *  input,
DATA_OBJ *  output,
double *  max 
)

Computes the root-sum-of-squares image from the input data.

input data can be either 3d: x-y-c, or 4d: x-y-z-c

The scaling factor needed to normalize the image (set the max pixel value to '1') is returned in 'max'

Parameters
inputinput data
outputoutput data (magnitude only, no phase)
maximage normalization scale factor (calculated, but not applied)
int compute_virtual_body_coil ( DATA_OBJ *  input,
DATA_OBJ *  output,
DATA_OBJ *  vbc_obj 
)

An implementation of [J], which builds on ideas in [K] for identifying alternate subspaces for reciever coil representations. This implementation normalizes each of the coils, and finds the dominant subspace support along the coil dimension of the data. This generates a virtual body coil, which can be used in parallel imaging to produce image reconstructions with more homogeneous illumination in a self-referenced way.

Usage:

With SENSE or SPACE RIP, one can renormalize the sensitivity estimates before the reconstruction step:

self_ref_b1_via_geyser( K, hdr, PhaseLines, W, mngr, 0 );
compute_virtual_body_coil( W, hdr.Nx, hdr.Ny, hdr.ncoils, Bv );
for (cnt=0;cnt<hdr.ncoils;cnt++) {
  for (cnt2=0;cnt2<hdr.Nx*hdr.Ny;cnt2++) {
     W[ cnt*hdr.Nx*hdr.Ny + cnt2 ] =
        complex_division( W[ cnt*hdr.Nx*hdr.Ny + cnt2 ],
                          Bv[ cnt2 ] );
  }
}

Alternatively, after reconstructing with GRAPPA, one can find the associated virtual body coil image via the projection across all coil images. See the grappa.c MEX file for an example.

Note that the virtual body coil image contains phase information, so it can be used in applications were a root-sum-of-squares body coil estimate is unsuitable.

Parameters
inputthe input coil data, sized [x-y-c] or [x-y-z-c]
outputthe output virtual body coil. If (output->data == NULL), memory will be allocated.
vbc_objreturns the coil combination coeficients, if non-NULL
int export_grappa_coef ( GRAPPA_PARMS **  gparm,
DATA_OBJ *  gexport,
int  ncoil,
int  nslc 
)

Converts an array of grappa parameters into a DATA_OBJ suitable for output to a .nd file, using write_nd_data2. Each location in the array is assumed to be for a specific slice of data, and the grappa kernel is presumed to be the same size (resulting in the same number of grappa coefficients) for each slice.

The returned data object will have the meta-data field populated, with a brief description of how the coefficients are arranged. This meta-data can be decoded in matlab using load_gparms.m or in C using import_grappa_coef(), with the data properly organized so that previously generated coefficients can be correctly used in future applications.

Parameters
gparmthe input array of grappa parameters, one set for each slice
gexportthe output data object
ncoilthe number of coils associated with the grappa parameters
nslcthe number of slices in the input array
void generate_b1_filter ( int  fill_fac,
short  Ny,
short  Nz_proc,
float *  samp_dens,
float *  filt_b1 
)

Generate a 2D Gaussian filter to extract the ky-kz region, near k-space center, to be used for sensitivity mapping.

Parameters
samp_densSampling density
filt_b1Gaussian filter to select ky-kz region to be used for B1 mapping.
int grappa ( COMPLEX *  I,
int *  PhaseLines,
dataheader *  hdr,
COMPLEX *  O,
int *  Nlines,
int **  OutPhaseLines,
int  YKS,
int  XKS,
float  radius,
int  ndk,
int *  dks_in 
)

An implementation of ``Generalized autocalibrating partially parallel acquisitions,'' by M. A. Griswold, P. M. Jakob, et. al.. Mag Reson Med, 47(6):1202-1210. Jun 2002. [DOI]

Processing is done in k-space, coil-by-coil. The reconstructed data includes all coils. To form an image, one can use root-sum-of-squares.

Parameters
Iinput data
PhaseLinesinput phase index list
hdrheader associated with input data
Obuffer to hold output data. This should be Nvertical x xres x Ncoils, i.e. zero padded
Nlinesnumber of non-zero output lines
OutPhaseLineslist of output non-zero lines (will be malloced)
YKSgrappa-kernel size, y-dim (2 or 4)
XKSgrappa-kernel size, x-dim (must be odd)
radiussize of reconstruction radius along freq-encode dimension. In some cases (e.g. coil sensitivity estimation) the full read-out length is not needed
ndknumber of Delta k's to compute. If zero, then the dks list is computed automatically.
dks_inthe list of Delta k's to compute. Set to NULL if ndk = 0.
int grappa_calc_parms ( dataheader *  hdr,
COMPLEX *  O,
int  YKS,
int *  indx0,
int  XKS,
int *  indx3,
int  nacs,
int *  acsl,
int  kxnum,
COMPLEX **  gparm,
int  fullsqr,
int  direction 
)

Calculate GRAPPA reconstruction parameters for the given data.

called by grappa

Example: for a 2x5 GRAPPA kernel, to reconstruct data sampled at 2x acceleration... YKS = 2; indx0 = {-1,1}; XKS = 5; indx3 = {-2,-1,0,1,2};

Parameters
hdrheader associated with input data
Ozero-padded k-space data
YKSsize of kernel along phase encode direction (length of indx0)
indx0phase encode offsets for kernel
XKSsize of kernel along read-out direction (length of indx3)
indx3sample offsets along readout
nacsnumber of ACS lines
acslarray of ACS line indices
kxnumnumber of readout (k_x) points used for calibration
gparmthe output GRAPPA parameters. Pass a NULL pointer, memory for the parameters will be allocated and returned.
fullsqrfor GROG: set to 1 if full dx set is needed, 2 for both dx and dy data set. set to 0 otherwise
directionfor GROG: modifies the yoffset value by {0 [default],+1,-1}.
int grappa_calc_recon_coef ( COMPLEX *  O,
int *  PhaseLines,
dataheader *  hdr,
GRAPPA_PARMS *  gparm,
float  radius 
)

The primary function for calculating GRAPPA coefficients.

Input the source calibration data (ACS) and the phase lines where that data was acquired, and the parameters will be returned in 'gparm'.

Parameters
Obuffer holding zero-padded input data. This should be Nvertical x xres x Ncoils
PhaseLinesinput list of phase index showing which lines of O have data.
hdrheader associated with input data
gparmneed to predeclare the kernel size (gparm->YKS and gparm->XKS) and gparm->parray (set to NULL). ndk and dks may also be preset, or set to ndk=0 and dks=NULL for values computed from PhaseLines. For standard GRAPPA parameters, set hdr->grog = 0; For GROG, use hdr->grog = 1 for a square parameter set, suitable for GROG along the phase encode direction. use hdr->grog = 2 for an additional set of parameters suitable for GROG along the readout direction.
radiuswidth, (0.0 .. 1.0], along the readout line to use in calibration
int grappa_for_cine ( short  Nx,
short  N_acq,
short  Ny,
int  Ny_proc,
short  Nz_proc,
short  Nf,
short  ncoils,
short  fill_fac,
short  fr1,
int  first_line,
float  pF,
KCOORD *  kline_c,
float *  samp_dens,
float *  weight_pF,
COMPLEX *  buf2,
COMPLEX *  b1,
COMPLEX *  data,
COMPLEX *  full 
)

A function to call GRAPPA from cine_unf

Parameters
NxSize of input dataset, along kx (readout)
N_acqk locs acquired / frame (in ky-kz plane)
NySize of input dataset, along ky (phase)
Ny_procNy used for processing, in vdsense
Nz_procNz used for processing
Nfnumber of frames (either time or phase)
ncoilsNumber of coils
fill_facfill factor, to accomodate non-int samp pts
fr1whether 1st frame considered even (0) or odd (1)
first_lineFirst y line kept, when cropping recon FOV
pFpartial-Fourier (1=fully sampled, 0.5=50% sampled)
kline_carray of acquired phase encode lines
b1pointer to coil sensitivity map data
datapointer to input data
fullpointer to reconstructed data
int grappa_recon ( dataheader *  hdr,
COMPLEX *  O,
int  YKS,
int *  indx0,
int  XKS,
int *  indx3,
int *  bins,
COMPLEX *  gparm,
float  radius 
)

reconstruct the missing data using GRAPPA.

called by grappa

Parameters
hdrstructure describing the data
Othe (zero-padded) input data
YKSsize of kernel along phase encode direction
indx0phase encode offsets for kernel
XKSsize of kernel along read-out direction
indx3sample offsets along readout
binsan array dscribing the data along the phase encoding dimension. 0 : empty ; 1 : raw data ; 2 : ACS line
gparmpointer to the GRAPPA reconstruction parameters
radiussize of reconstruction radius along freq-encode dimension
int grog_generate_eigd ( GRAPPA_PARMS **  gparm,
int  ncoil,
int  nslc 
)

Derives the eigen-vectors and eigen-values needed by GROG from the input GRAPPA parameters.

These are stored in the .eig_d and .eig_v pointers of the GRAPPA_OBJ, over-writing whatever may have been there previously.

Parameters
gparminput array of GRAPPA parameters, suitable for GROG ( set gparm->grog = 1 or 2 prior to generating gparm coefficients )
ncoilnumber of coils assocated with the input array
nslcnumber of slices assocated with the input array
int grog_grid_2Ddata ( DATA_OBJ *  smpl_x,
DATA_OBJ *  smpl_y,
DATA_OBJ *  kdata,
DATA_OBJ *  src,
DATA_OBJ *  new_x,
DATA_OBJ *  new_y,
DATA_OBJ *  kdata_out,
GRAPPA_PARMS *  gparms 
)

A function to shift k-space data locations, using GROG.

the input data is a triplet array, [ smpl_x, smpl_y, kdata ] that is transformed to [ new_x, new_y, kdata_out ]. The first 5 elements are specified, while the 6th is computed.

Parameters
smpl_xa 1D (double / "dblr") array of sample point locations, along the x (horizontal) axis
smpl_ya 1D (double / "dblr") array of sample point locations, along the y (vertical) axis
kdataa 1D (complex double / "dblc") array of sampled k-space data, at each {smpl_x,smpl_y} point in the previous two arrays
srca 1D (short int / "intr") array of indices, indicating which kdata point should be used for each output target listed in {new_x,new_y}
new_xa 1D (double / "dblr") array of destination grid locations, along x
new_ya 1D (double / "dblr") array of destination grid locations, along y
kdata_outa 1D (complex double / "dblc") array of computed k-space data, at the points specified by {new_x,new_y}
gparmsa pointer to the GROG parameters
int import_grappa_coef ( GRAPPA_PARMS ***  gparm_out,
DATA_OBJ *  gimport 
)

Converts a DATA_OBJECT, with an approriately defined and consistent metadata field, into an array of grappa parameters suitable for subsequent GRAPPA reconstruction calls. Each location in the array is assumed to be for a specific slice of data, and the grappa kernel is presumed to be the same size (resulting in the same number of grappa coefficients) for each slice.

Parameters
gparm_outa pointer to an (empty) array of GRAPPA parameters. This function will replace the data pointer with newly allocated memory.
gimportthe input data object
int inv_sens_matrix ( short  ncoils,
short  npts,
float  thresh,
double **  a,
double *  w,
double **  v,
double *  wgt_roe,
double *  wgt,
double **  S,
double **  S_buf,
double **  S_inv,
double **  S_inv_buf 
)

Estimate the inverse of the sensitivity matrix, for each pixel.

The threshold paramter is described in Madore B. "UNFOLD-SENSE: A parallel MRI method with self-calibration and artifact suppression." Magn Reson Med 2004; 52(2):310–320.

At the expense of artifact content (presumed to be suppressed with UNFOLD later in the processing), a matrix with less noise (i.e. greater regularization) for the ith entry will be obtained. If the npts vectors representing how each one of the overlapped points is seen by the coils were orthogonal, optimum SNR (i.e. the Roemer case) would be achieved. Since the solution obtained by matrix inversion would amplify noise beyond the threshold, fool the recon to believe these vectors are more orthogonal than they really are, so that it finds smaller weights. The vector for the ith component, the one under consideration, is left unchanged, while the npts-1 others are made more orthogonal to it. Each one of these npts-1 vectors has a component parallel, and one orthogonal to the ith (unit) vector uv[]. Vectors are made more orthogonal to uv[] by keeping only a fraction of the parallel component (and leaving the orthogonal one untouched). The goal of the iterative process is to find the right value that brings noise amplification down to a factor of thresh.

The implication is that the reconstructed image will have less noise than a standard psuedo-inverse, at the expense of greater artifact. It is assumed that these artifacts will then be eliminated using UNFOLD.

If UNFOLD is not being used, then one can use a very large value for thresh to suppress these iterations.

Parameters
ncoilsNumber of receiver coils.
nptsNumber of non-degenerate points (npts <= n_s, the acceleration seen by SENSE).
threshNoise threshold.
aMatrix to invert for the SENSE recon of npts pixels. (of size 2*ncoils x 2*accel
wSingular values, in the SENSE recon. (of size 2*accel)
vRequired for backsubstitution, in the SENSE recon. (of size (2*accel)
wgt_roeWeights as in Roemer et al.
wgtWeights with present method.
SSensitivity matrix, copy of "a".
S_bufModified version of S.
S_invInverse of the sensitivity matrix. (2*accel x 2*ncoils)
S_inv_bufUnfinished version of S_inv.
int rlsgrappa ( COMPLEX *  I,
int *  PhaseLines,
dataheader *  hdr,
COMPLEX *  O,
int *  Nlines,
int **  OutPhaseLines,
int  YKS,
float  radius,
int  ndk,
int *  dks_in,
float  lambda 
)

An implementation of ``RLS-GRAPPA: Reconstructing parallel MRI data with adaptive filters,'' by W S Hoge, F Gallego, Z Xiao, and D H Brooks, to appear in Proc 2008 IEEE Intl Symp on Biomedical Imaging, May 2008.

Processing is done in hybrid-space, ky-x, coil-by-coil with GRAPPA reconstruction coefficients calculated simlutaneously with missing data reconstruction using an adapative RLS filter. Currently only 'RLS-along-x' is supported.

Parameters
Iinput data
PhaseLinesinput phase index list
hdrheader associated with input data
Obuffer to hold output data. This should be Nvertical x xres x Ncoils, i.e. zero padded
Nlinesnumber of non-zero output lines
OutPhaseLineslist of output non-zero lines (will be malloced)
YKSgrappa-kernel size, y-dim (2 or 4)
radiussize of reconstruction radius along freq-encode dimension. In some cases (e.g. coil sensitivity estimation) the full read-out length is not needed
ndknumber of Delta k's to compute. If zero, then the dks list is computed automatically.
dks_inthe list of Delta k's to compute. Set to NULL if ndk = 0.
lambdathe rls forgetting factor
void self_ref_b1_basic ( COMPLEX *  kspace_data,
float  b1snr_thr,
short  Ny,
short  Nx,
short  Nz,
short  ncoils,
short  Nf,
fil_fft_mngr *  mngr,
COMPLEX *  b1 
)

Calculate B1 sensitivity maps from raw kspace data.

This is a very basic function: it applies only a 1/2 sinusoidal filter only along ky; assumes kspace is fully sampled with the DC coordinate in the center of k-space. An fftshift is during processing to ensure that the estimated coil sensitivities have smooth phase.

Parameters
kspace_datakspace data of size (ky, kx, z, coils, t).
b1snr_thrSNR threshold for sensitivity maps.
NyDimension along the y/ky direction.
NxDimension along the x/kx direction.
NzDimension along the z/kz direction.
ncoilsNumber of coil elements in array.
NfNumber of frames in kspace_data.
mngrFFT manager
b1B1 sensitivity maps. Size: nx-ny-nz-ncoils
void self_ref_b1_for_cine ( COMPLEX *  kspace_data,
float  sigmak,
float  b1snr_thr,
short  N_acq,
short  Nx,
short  Ny,
short  Ny_proc,
short  Nz_proc,
short  ncoils,
short  Nf,
short  fr1,
short  fill_fac,
short  first_line,
float *  filt_b1,
KCOORD *  kline_o,
KCOORD *  kline_e,
float *  samp_dens,
fil_fft_mngr *  mngr,
COMPLEX *  bufsp,
COMPLEX *  b1 
)

Calculate B1 sensitivity maps from the dynamic data itself (i.e. self- referenced). It uses a combination of the strategies as part of GRAPPA (k- space center sampled more densely) and TSENSE (UNFOLD-based approach).

Parameters
kspace_dataRe-ordered data of size (Nx, N_acq, Nz, nph, ncoils).
sigmakNoise level in kspace_data, as determined by eval_knoise_cine.
b1snr_thrSNR threshold for sensitivity maps.
N_acqNumber of sampled lines per frame.
NxDimension along the x/kx direction.
NyDimension along the y/ky direction.
Ny_proc

Ny if sampling density > 1 at ky~=0.

Nz_procDimension along the z/kz direction.
ncoilsNumber of coil elements in array.
NfNumber of frames in kspace_data.
fr11st frame considered even (0) or odd (1)
fill_facFactor to make ky lines integer.
first_lineFirst line for cropping (0 = off).
filt_b1Filter to select ky-kz region for B1 sensitivity mapping.
kline_oSampling function, odd frames.
kline_eSampling function, even frames.
samp_densSampling density, in ky-kz plane.
mngrFFT manager
bufspMemory space, to work in.
b1calculated B1 sensitivity maps.
int self_ref_b1_via_geyser ( COMPLEX *  datain,
dataheader *  hdr,
int *  PhaseLines,
COMPLEX *  ptrW,
fil_fft_mngr *  mngr,
int  m 
)

Compute the coil sensitivity estimate for one frame of data, using GRAPPA as an intermediary to fill in any gaps (of $ 2 \Delta k $) in the low-frequency coordinate range of the acquired k-space data.

For reference, see Hoge WS and Brooks DH, "Using GRAPPA to improve auto-calibrated coil sensitivity estimation for the SENSE family of parallel imaging reconstruction algorithms," Magn. Reson. Med., Magn Reson Med, 2008; 60(2):462-467.

Parameters
dataininput k-space data (size: P-N-L)
hdra header to describe the data and associated output matrix size
PhaseLinesa list of acquired phase encode indices, of length P
ptrWA buffer to store output coil sensitivity maps. (size: M-N-L)
mngra pointer to an FFT manager
mbinary flag: use a mask? default: 0
int spacerip_for_cine ( short  Nx,
short  N_o,
short  N_e,
short  N_c,
short  Ny,
int  Ny_proc,
short  Nz_proc,
short  Nf,
short  ncoils,
float  pF,
KCOORD *  kline_o,
KCOORD *  kline_e,
KCOORD *  kline_c,
COMPLEX *  b1,
COMPLEX *  buf,
COMPLEX *  data,
COMPLEX *  spacerip,
short  recon_flag 
)

An implementation of "Sensitivity profiles from an array of coils for encoding and reconstruction in parallel (SPACE RIP)." W. E. Kyriakos, et. al. Magn Reson Med. 44(2):301-308, August 2000.

SPACE RIP can be used to reconstruct data sampled non-uniformly on a rectilinear sampling grid.

Parameters
NxSize of input dataset, along kx (readout)
N_ok locs acquired in odd ky-kz frames
N_ek locs acquired in even ky-kz frames
N_ck locs acquired in combined frames
NyFull matrix size along ky (phase)
Ny_procNy used for processing, in vdsense
Nz_procNz used for processing
Nfnumber of frames (either time or phase)
ncoilsNumber of coils
pFpartial-Fourier (1=fully sampled, 0.5=50% sampled)
kline_olines acquired on 'odd' frames
kline_elines acquired on 'even' frames
kline_clines acquired on combined frames
b1pointer to sensitivity map data
bufpointer to available memory block
datapointer to input data
spacerippointer to reconstructed data
recon_flag0 for low t freqs, 1 for high t freqs
int sprip3d_conv_pts2y ( int  nph,
int *  phz_y,
int *  phz_z,
int  M,
int  N,
int *  Y 
)

calculate the array index values for the specified acquired phase encode points. i.e. maps (phz_y,phz_z) = (2,1) of a (MxN = ) 10 x 8 image to y = 12

y[i] = phz_z[i] * M + phz_y[i]

  01234567
0 ........
1 ........
2 .x......
3 ........
4 ........
5 ........
6 ........
7 ........
8 ........
9 ........

Note: the C convention (which is followed here) is to start counting at zero, which differs from Matlab, which starts counting vector elements at one.

Parameters
nphnumber of acquired phase encode points
phz_yy locs of each phase encode point
phz_zz locs of each phase encode point
Mtotal span of the FOV, along y
Ntotal span of the FOV, along z
Ycalcuated list of array index values
void sprip3d_lsqr_init ( int  nph,
int  M,
int  N,
int  L 
)

setup the memory allocation and FFT configuration for the LSQR reconstruction algorithm of 3D imaging data, sprip3d_lsqrsolv().

Parameters
nphnumber of acquired phase encode lines
Mtotal span of the phase-encode dimension
Ntotal span of the phase-encode dimension
Lnumber of coils
void sprip3d_lsqr_post ( )

tear-down the memory allocation and FFT configuration used during the LSQR reconstructions of 3D data.

void sprip3d_lsqrsolv ( int *  Y,
int  nph,
int  m,
int  n,
int  L,
COMPLEX *  W,
COMPLEX *  Sk,
COMPLEX **  I,
float *  inlams,
int  lamsize,
int  maxits,
double *  xnormout,
double *  Rout,
int  tid 
)

Solves the inverse problem associated with one slice of the 3D SPACE RIP reconstruction algorithm, utilizing fast matrix vector products available via FFT operations and embedded Krylov techniques to quickly find solutions associated with multiple regularization values.

Parameters
Yphase encode index list
nphnumber of phase encodes
msize of full phase encode span along y
nsize of full phase encode span along z
Lnumber of coils
W(M,N,L) coil sensitivity data
Sk(nph,L) acquired data
Ipointer to reconstructed image data. If this value points to an existing memory block, a single image is returned, corresponding to the highest regularization parameter value with no aliasing in the image. If a NULL pointer is passed, then all reconstructions (one for each lambda) is returned. After copying the most desirable image, free the buffer to prevent memory leaks.
inlamslist of lambda's
lamsizenumber of lambda's in the lams list
maxitsmax number of iterations
xnormoutreturns the norm of solution, for each lambda
Routreturns the norm of residual error, for each lambda
tidunused
void sprip_cgsolv ( int *  Y,
int  nph,
int  n,
int  L,
COMPLEX *  Wx0,
COMPLEX *  u,
COMPLEX *  x,
COMPLEX  lambda,
int  k_max,
int  tid 
)

Solves the inverse problem associated with one column of the SPACE RIP reconstruction algorithm, utilizing fast matrix vector products available via FFT operations

Parameters
Yphase encode index list
nphnumber of phase encodes
nsize of full phase encode span
Lnumber of coils
Wx0(n,L) coil sensitivity data
u(nph,L) acquired data
x(length n) reconstructed image data
lambdaDLS regularization setting
k_maxmax number of iterations
tidthread id
void sprip_cgsr_init ( int  nph,
int  N,
int  L,
int  threadnum 
)

Set up memory allocation blocks and FFT plans needed for cgsolv

Parameters
nphnumber of acquired phase encode lines
Ntotal span of the phase-encode dimension
Lnumber of coils
threadnumnumber of threads to prepare for
void sprip_cgsr_post ( )

Tear-down the allocated memory blocks and FFT mananger used with cgsolv

void sprip_cgsr_recon ( float *  phz_encode_list,
int  P,
int  N,
int  M,
int  L,
COMPLEX *  W,
COMPLEX *  Sk,
double  tau,
int  maxits,
COMPLEX *  rho,
int  tmax 
)

Reconstruct an image using the Space-RIP paradigm, solving the linear system via the method of conjugate gradients

Parameters
phz_encode_listthe phase encodes used (length P vector) range = [0, N-1], < with the center of k-space at N/2
Pnumber of phase encodes used
Nnumber of reconstructed samples along phase encode direction
Mnumber of samples along readout direction
Lnumber of coils
Wcoil sensitivity estimate in spatial domain of size (N,M,L)
Skacquired data in k-space domain of size (P,M,L)
tauthe regularization parameter
maxitsmax number of iterations
rhothe reconstructed image (size = M-by-N)
tmaxmax number of threads to use. valid range: [0..M] If 0, then default number is used. set to non-zero if memory workspace is pre-allocated.
void sprip_lsqr_init ( int  nph,
int  N,
int  L,
int  threadnum,
int  lamsize,
int  maxits 
)

setup the memory allocation and FFT manager for the LSQR reconstruction algorithm, lsqrsolv.

Parameters
nphnumber of acquired phase encode lines
Ntotal span of the phase-encode dimension
Lnumber of coils
threadnumnumber of threads to prepare for
lamsizemaximum number of lambda's
maxitsmaximum number of iterations
void sprip_lsqr_post ( )

tear-down the memory allocation and FFT manager used during the LSQR reconstructions.

void sprip_lsqr_recon ( float *  phz_encode_list,
int  P,
int  N,
int  M,
int  L,
COMPLEX *  W,
COMPLEX *  Sk,
float *  lams,
int  numlams,
COMPLEX *  D,
int  maxits,
COMPLEX **  rho,
double *  xout,
double *  rout,
int  tmax,
float *  usedlams 
)

reconstruct an image using the SPACE-RIP paradigm, solving the linear system using the LSQR-Hybrid algorithm. This algorithm embeds the regularization parameter search within each iteration, efficiently producing multiple solutions for each regularization paramater (lambda) value.

Efficiently solves:

\[ \min_x \{ \| A x - b \|_2 + \lambda \| D x \|_2 \} \]

for multiple vales of $ \lambda $, where $ D $ is a diagonal regularization operator matrix.

Parameters
phz_encode_listthe phase encodes used (length P vector) range = [0, N-1], with the center of k-space at N/2
Pnumber of phase encodes used
Nnumber of reconstructed samples along phase encode direction
Mnumber of samples along readout direction
Lnumber of coils
Wcoil sensitivity estimate in spatial domain of size ( N, M, L )
Skacquired data in k-space domain of size ( P, M, L, nimgs )
lamsarray of regularization parameter. values (if *lams == NULL, use the default of 10^{-4} to 10^1 on a logarithmic scale.)
numlamsnumber of lambda values in the array (max=50). if zero, then default=20 is used
Da length N vector for varied regularization along the reconstructed dimension. Set to NULL to use the identity matrix.
maxitsthe maximum number of iterations
rhothe reconstructed image(s). If given a pointer, then the algorithm will return one image of size M-by-N. if (*rho==NULL), it will return an image for each of the regularization values. To avoid memory leaks, clear the returned buffer after the most desirable image is copied.
xoutreturns the estimated norm of the solution for each lambda
routreturns the norm of residual error, for each lambda
tmaxmaximum number of threads to use
usedlamsan array of size M, used to return the regularization value determined by the LSQR-Hybrid algorithm for each solved problem.
void sprip_lsqrsolv ( int *  Y,
int  nph,
int  N,
int  L,
COMPLEX *  Wx0,
COMPLEX *  Sk,
COMPLEX **  I,
float *  lams,
int  lamsize,
COMPLEX *  D,
int  maxits,
double *  xnormout,
double *  Rout,
int  tid,
float *  lout 
)

Solves the inverse problem associated with one column of the SPACE RIP reconstruction algorithm, utilizing fast matrix vector products available via FFT operations and embedded Krylov techniques to quickly find solutions associated with multiple regularization values.

Parameters
Yphase encode index list
nphnumber of phase encodes
Nsize of full phase encode span
Lnumber of coils
Wx0(N,L) coil sensitivity data
Sk(nph,L) acquired data
Ireconstructed image data
lamslist of lambda's
lamsizenumber of lambda's in the lams list
Da vector to use in non-I regularization. Set to NULL if not needed/desired.
maxitsmax number of iterations
xnormoutreturns the norm of solution, for each lambda
Routreturns the norm of residual error, for each lambda
tidthread id
void vdsense ( short  Nx,
short  N_acq,
short  Ny,
short  Ny_proc,
short  Nz_proc,
short  Nf,
short  ncoils,
short  fill_fac,
short  fr1,
short  first_line,
float  pF,
fil_fft_mngr *  mngr,
KCOORD *  kline_o,
KCOORD *  kline_e,
float *  samp_dens,
float *  weight_pF,
COMPLEX *  buf,
DATATYPE *  Uarray,
COMPLEX *  data,
COMPLEX *  vdsen_ima 
)

The reconstruction strategy implemented by Variable density SENSE (VD-SENSE) is to apply a Cartesian SENSE reconstruction to data that is not necessarily uniformly sub-sampled. The processing proceeds in three steps.

First, the coil sensitivity estimates are used to construct reconstruction (or un-mixing or un-aliasing) operators, based on the highest acceleration factor employed in the sampling scheme. Second, the raw data is zero padded and transformed to the spatial domain. Finally, the reconstruction operators are applied to the spatial domain representation of the acquired image data, to form the un-aliased image.

Parameters
NxSize in x direction
N_acqnumber of acquired samples along y
NySize in y direction
Ny_procNy used for processing
Nz_procNz used for processing (number of slices)
NfSize in time direction
ncoilsnumber of coils
fill_facfill factor, to accomodate non-int samp pts
fr1whether 1st frame considered even (0) or odd (1)
first_lineFirst y line kept, when cropping recon FOV
pFpartial-Fourier (1=fully sampled, 0.5=50% sampled)
mngra pointer to an FFT manager
UarrayContains an inverted SENSE matrix U for each pixel, to perform the unaliasing.
void vdsense_apply ( short  Nx,
short  Ny,
short  Nz,
short  Nf,
short  ncoils,
COMPLEX *  sen_before,
DATATYPE *  Uarray,
COMPLEX *  sen_after 
)

Apply variable-density SENSE.

Parameters
NxSize in x direction.
NySize in y direction.
NzSize in z direction.
NfSize time direction.
ncoilsNumber of receiver coils.
sen_beforeData, before vdsense. (size: Nx*Ny*Nz*Nf*ncoils*sizeof(COMPLEX) )
UarrayContains an inverted matrix U for each pixel, to unalias.
sen_afterData, after vdsense. (size: Nx*Ny*Nz*Nf*sizeof(COMPLEX))
void vdsense_prep ( short  Nx,
short  Ny,
short  Nz,
short  ncoils,
float  thresh,
int  accel,
XCOORD *  alias,
COMPLEX *  b1,
DATATYPE *  Uarray 
)

Invert sensitivity matrices, in preparation for a vdsense reconstruction.

See inv_sens_matrix() for a description of the thresh parameter.

Parameters
NxSize in x direction.
NySize in y direction.
NzSize in z direction.
ncoilsNumber of receiver coils.
threshThreshold for noise tolerance.
accelthe maximum delta k in the sampling pattern
aliasDescription of the aliasing PSF (where aliasing comes from).
b1B1 sensitivity maps. (size: Nx*Ny*Nz*sizeof(COMPLEX))
UarrayArray where all the inverted matrices will be stored. (size: Nx*Ny*Nz*ncoils*sizeof(COMPLEX))