ncigt-fil
|
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) |
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]:
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.
where 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.
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
Pruessmann KP, Weiger M, Scheidegger MB, Boesiger P. SENSE: Sensitivity encoding for fast MRI. Magn Reson Med 1999; 42(5):952-62.
Madore B. Using UNFOLD to remove artifacts in parallel imaging and in partial-Fourier imaging. Magn Reson Med 2002; 48(3):493-501.
Kellman P, Epstein FH, McVeigh ER. Adaptive sensitivity encoding incorporating temporal filtering (TSENSE). Magn Reson Med 2001; 45(5):846-852.
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.
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.
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.
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.
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
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.
Buehrer M, Boesiger P, and Kozerke S. "Virtual body coil calibration for phased-array imaging." 17th ISMRM Scientific Meeting, page 759, Honolulu, HI, 2009.
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.
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
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.
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.
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
input | the input coil data, sized [x-y-c] or [x-y-z-c] |
vbc_obj | returns 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'
input | input data |
output | output data (magnitude only, no phase) |
max | image 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.
input | the input coil data, sized [x-y-c] or [x-y-z-c] |
output | the output virtual body coil. If (output->data == NULL), memory will be allocated. |
vbc_obj | returns 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.
gparm | the input array of grappa parameters, one set for each slice |
gexport | the output data object |
ncoil | the number of coils associated with the grappa parameters |
nslc | the 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.
samp_dens | Sampling density |
filt_b1 | Gaussian 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.
I | input data |
PhaseLines | input phase index list |
hdr | header associated with input data |
O | buffer to hold output data. This should be Nvertical x xres x Ncoils, i.e. zero padded |
Nlines | number of non-zero output lines |
OutPhaseLines | list of output non-zero lines (will be malloced) |
YKS | grappa-kernel size, y-dim (2 or 4) |
XKS | grappa-kernel size, x-dim (must be odd) |
radius | size of reconstruction radius along freq-encode dimension. In some cases (e.g. coil sensitivity estimation) the full read-out length is not needed |
ndk | number of Delta k's to compute. If zero, then the dks list is computed automatically. |
dks_in | the 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};
hdr | header associated with input data |
O | zero-padded k-space data |
YKS | size of kernel along phase encode direction (length of indx0) |
indx0 | phase encode offsets for kernel |
XKS | size of kernel along read-out direction (length of indx3) |
indx3 | sample offsets along readout |
nacs | number of ACS lines |
acsl | array of ACS line indices |
kxnum | number of readout (k_x) points used for calibration |
gparm | the output GRAPPA parameters. Pass a NULL pointer, memory for the parameters will be allocated and returned. |
fullsqr | for GROG: set to 1 if full dx set is needed, 2 for both dx and dy data set. set to 0 otherwise |
direction | for 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'.
O | buffer holding zero-padded input data. This should be Nvertical x xres x Ncoils |
PhaseLines | input list of phase index showing which lines of O have data. |
hdr | header associated with input data |
gparm | need 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. |
radius | width, (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
Nx | Size of input dataset, along kx (readout) |
N_acq | k locs acquired / frame (in ky-kz plane) |
Ny | Size of input dataset, along ky (phase) |
Ny_proc | Ny used for processing, in vdsense |
Nz_proc | Nz used for processing |
Nf | number of frames (either time or phase) |
ncoils | Number of coils |
fill_fac | fill factor, to accomodate non-int samp pts |
fr1 | whether 1st frame considered even (0) or odd (1) |
first_line | First y line kept, when cropping recon FOV |
pF | partial-Fourier (1=fully sampled, 0.5=50% sampled) |
kline_c | array of acquired phase encode lines |
b1 | pointer to coil sensitivity map data |
data | pointer to input data |
full | pointer 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
hdr | structure describing the data |
O | the (zero-padded) input data |
YKS | size of kernel along phase encode direction |
indx0 | phase encode offsets for kernel |
XKS | size of kernel along read-out direction |
indx3 | sample offsets along readout |
bins | an array dscribing the data along the phase encoding dimension. 0 : empty ; 1 : raw data ; 2 : ACS line |
gparm | pointer to the GRAPPA reconstruction parameters |
radius | size 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.
gparm | input array of GRAPPA parameters, suitable for GROG ( set gparm->grog = 1 or 2 prior to generating gparm coefficients ) |
ncoil | number of coils assocated with the input array |
nslc | number 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.
smpl_x | a 1D (double / "dblr") array of sample point locations, along the x (horizontal) axis |
smpl_y | a 1D (double / "dblr") array of sample point locations, along the y (vertical) axis |
kdata | a 1D (complex double / "dblc") array of sampled k-space data, at each {smpl_x,smpl_y} point in the previous two arrays |
src | a 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_x | a 1D (double / "dblr") array of destination grid locations, along x |
new_y | a 1D (double / "dblr") array of destination grid locations, along y |
kdata_out | a 1D (complex double / "dblc") array of computed k-space data, at the points specified by {new_x,new_y} |
gparms | a 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.
gparm_out | a pointer to an (empty) array of GRAPPA parameters. This function will replace the data pointer with newly allocated memory. |
gimport | the 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.
ncoils | Number of receiver coils. |
npts | Number of non-degenerate points (npts <= n_s, the acceleration seen by SENSE). |
thresh | Noise threshold. |
a | Matrix to invert for the SENSE recon of npts pixels. (of size 2*ncoils x 2*accel |
w | Singular values, in the SENSE recon. (of size 2*accel) |
v | Required for backsubstitution, in the SENSE recon. (of size (2*accel) |
wgt_roe | Weights as in Roemer et al. |
wgt | Weights with present method. |
S | Sensitivity matrix, copy of "a". |
S_buf | Modified version of S. |
S_inv | Inverse of the sensitivity matrix. (2*accel x 2*ncoils) |
S_inv_buf | Unfinished 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.
I | input data |
PhaseLines | input phase index list |
hdr | header associated with input data |
O | buffer to hold output data. This should be Nvertical x xres x Ncoils, i.e. zero padded |
Nlines | number of non-zero output lines |
OutPhaseLines | list of output non-zero lines (will be malloced) |
YKS | grappa-kernel size, y-dim (2 or 4) |
radius | size of reconstruction radius along freq-encode dimension. In some cases (e.g. coil sensitivity estimation) the full read-out length is not needed |
ndk | number of Delta k's to compute. If zero, then the dks list is computed automatically. |
dks_in | the list of Delta k's to compute. Set to NULL if ndk = 0. |
lambda | the 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.
kspace_data | kspace data of size (ky, kx, z, coils, t). |
b1snr_thr | SNR threshold for sensitivity maps. |
Ny | Dimension along the y/ky direction. |
Nx | Dimension along the x/kx direction. |
Nz | Dimension along the z/kz direction. |
ncoils | Number of coil elements in array. |
Nf | Number of frames in kspace_data. |
mngr | FFT manager |
b1 | B1 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).
kspace_data | Re-ordered data of size (Nx, N_acq, Nz, nph, ncoils). |
sigmak | Noise level in kspace_data, as determined by eval_knoise_cine. |
b1snr_thr | SNR threshold for sensitivity maps. |
N_acq | Number of sampled lines per frame. |
Nx | Dimension along the x/kx direction. |
Ny | Dimension along the y/ky direction. |
Ny_proc |
|
Nz_proc | Dimension along the z/kz direction. |
ncoils | Number of coil elements in array. |
Nf | Number of frames in kspace_data. |
fr1 | 1st frame considered even (0) or odd (1) |
fill_fac | Factor to make ky lines integer. |
first_line | First line for cropping (0 = off). |
filt_b1 | Filter to select ky-kz region for B1 sensitivity mapping. |
kline_o | Sampling function, odd frames. |
kline_e | Sampling function, even frames. |
samp_dens | Sampling density, in ky-kz plane. |
mngr | FFT manager |
bufsp | Memory space, to work in. |
b1 | calculated 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 ) 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.
datain | input k-space data (size: P-N-L) |
hdr | a header to describe the data and associated output matrix size |
PhaseLines | a list of acquired phase encode indices, of length P |
ptrW | A buffer to store output coil sensitivity maps. (size: M-N-L) |
mngr | a pointer to an FFT manager |
m | binary 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.
Nx | Size of input dataset, along kx (readout) |
N_o | k locs acquired in odd ky-kz frames |
N_e | k locs acquired in even ky-kz frames |
N_c | k locs acquired in combined frames |
Ny | Full matrix size along ky (phase) |
Ny_proc | Ny used for processing, in vdsense |
Nz_proc | Nz used for processing |
Nf | number of frames (either time or phase) |
ncoils | Number of coils |
pF | partial-Fourier (1=fully sampled, 0.5=50% sampled) |
kline_o | lines acquired on 'odd' frames |
kline_e | lines acquired on 'even' frames |
kline_c | lines acquired on combined frames |
b1 | pointer to sensitivity map data |
buf | pointer to available memory block |
data | pointer to input data |
spacerip | pointer to reconstructed data |
recon_flag | 0 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.
nph | number of acquired phase encode points |
phz_y | y locs of each phase encode point |
phz_z | z locs of each phase encode point |
M | total span of the FOV, along y |
N | total span of the FOV, along z |
Y | calcuated 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().
nph | number of acquired phase encode lines |
M | total span of the phase-encode dimension |
N | total span of the phase-encode dimension |
L | number 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.
Y | phase encode index list |
nph | number of phase encodes |
m | size of full phase encode span along y |
n | size of full phase encode span along z |
L | number of coils |
W | (M,N,L) coil sensitivity data |
Sk | (nph,L) acquired data |
I | pointer 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. |
inlams | list of lambda's |
lamsize | number of lambda's in the lams list |
maxits | max number of iterations |
xnormout | returns the norm of solution, for each lambda |
Rout | returns the norm of residual error, for each lambda |
tid | unused |
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
Y | phase encode index list |
nph | number of phase encodes |
n | size of full phase encode span |
L | number of coils |
Wx0 | (n,L) coil sensitivity data |
u | (nph,L) acquired data |
x | (length n) reconstructed image data |
lambda | DLS regularization setting |
k_max | max number of iterations |
tid | thread id |
void sprip_cgsr_init | ( | int | nph, |
int | N, | ||
int | L, | ||
int | threadnum | ||
) |
Set up memory allocation blocks and FFT plans needed for cgsolv
nph | number of acquired phase encode lines |
N | total span of the phase-encode dimension |
L | number of coils |
threadnum | number 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
phz_encode_list | the phase encodes used (length P vector) range = [0, N-1], < with the center of k-space at N/2 |
P | number of phase encodes used |
N | number of reconstructed samples along phase encode direction |
M | number of samples along readout direction |
L | number of coils |
W | coil sensitivity estimate in spatial domain of size (N,M,L) |
Sk | acquired data in k-space domain of size (P,M,L) |
tau | the regularization parameter |
maxits | max number of iterations |
rho | the reconstructed image (size = M-by-N) |
tmax | max 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
.
nph | number of acquired phase encode lines |
N | total span of the phase-encode dimension |
L | number of coils |
threadnum | number of threads to prepare for |
lamsize | maximum number of lambda's |
maxits | maximum 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:
for multiple vales of , where
is a diagonal regularization operator matrix.
phz_encode_list | the phase encodes used (length P vector) range = [0, N-1], with the center of k-space at N/2 |
P | number of phase encodes used |
N | number of reconstructed samples along phase encode direction |
M | number of samples along readout direction |
L | number of coils |
W | coil sensitivity estimate in spatial domain of size ( N, M, L ) |
Sk | acquired data in k-space domain of size ( P, M, L, nimgs ) |
lams | array of regularization parameter. values (if *lams == NULL, use the default of 10^{-4} to 10^1 on a logarithmic scale.) |
numlams | number of lambda values in the array (max=50). if zero, then default=20 is used |
D | a length N vector for varied regularization along the reconstructed dimension. Set to NULL to use the identity matrix. |
maxits | the maximum number of iterations |
rho | the 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. |
xout | returns the estimated norm of the solution for each lambda |
rout | returns the norm of residual error, for each lambda |
tmax | maximum number of threads to use |
usedlams | an 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.
Y | phase encode index list |
nph | number of phase encodes |
N | size of full phase encode span |
L | number of coils |
Wx0 | (N,L) coil sensitivity data |
Sk | (nph,L) acquired data |
I | reconstructed image data |
lams | list of lambda's |
lamsize | number of lambda's in the lams list |
D | a vector to use in non-I regularization. Set to NULL if not needed/desired. |
maxits | max number of iterations |
xnormout | returns the norm of solution, for each lambda |
Rout | returns the norm of residual error, for each lambda |
tid | thread 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.
Nx | Size in x direction |
N_acq | number of acquired samples along y |
Ny | Size in y direction |
Ny_proc | Ny used for processing |
Nz_proc | Nz used for processing (number of slices) |
Nf | Size in time direction |
ncoils | number of coils |
fill_fac | fill factor, to accomodate non-int samp pts |
fr1 | whether 1st frame considered even (0) or odd (1) |
first_line | First y line kept, when cropping recon FOV |
pF | partial-Fourier (1=fully sampled, 0.5=50% sampled) |
mngr | a pointer to an FFT manager |
Uarray | Contains 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.
Nx | Size in x direction. |
Ny | Size in y direction. |
Nz | Size in z direction. |
Nf | Size time direction. |
ncoils | Number of receiver coils. |
sen_before | Data, before vdsense. (size: Nx*Ny*Nz*Nf*ncoils*sizeof(COMPLEX) ) |
Uarray | Contains an inverted matrix U for each pixel, to unalias. |
sen_after | Data, 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.
Nx | Size in x direction. |
Ny | Size in y direction. |
Nz | Size in z direction. |
ncoils | Number of receiver coils. |
thresh | Threshold for noise tolerance. |
accel | the maximum delta k in the sampling pattern |
alias | Description of the aliasing PSF (where aliasing comes from). |
b1 | B1 sensitivity maps. (size: Nx*Ny*Nz*sizeof(COMPLEX)) |
Uarray | Array where all the inverted matrices will be stored. (size: Nx*Ny*Nz*ncoils*sizeof(COMPLEX)) |