Comparison of MRI Under-Sampling Techniques for Compressed Sensing with Translation Invariant Wavelets Using FastTestCS: A Flexible Simulation Tool

Show more

1. Introduction

Magnetic Resonance Imaging (MRI) is a diagnostic modality used to create in-vivo images of 3-Dimensional (3D) biological tissue utilizing magnetic fields, gradients and receivers. When MR signal k-space data are fully sampled based on the Nyquist sampling criteria, some typical high resolution scans take five minutes, allowing time for patient and biologic movement, which negatively impacts Image Quality (IQ). Therefore, it is of great desire to accelerate the data acquisition and achieve better diagnostic images of the body. A promising theory is Compressed Sensing (CS) to under-sample k-space below what the Nyquist criteria requires without compromising IQ. Candès et al. introduce theory and experiments for the general CS problem [1] . Lustig et al. apply CS to the MRI setting with in-depth and practical design descriptions and experiments [2] .

Under-sampling for CS must be carefully designed so that artifacts and errors can be minimized and incoherent. Parallel imaging is another speed-up technique that combines the signals from multiple coils at the same time allowing some under-sampling of k-space to boost the IQ. Deshmane et al. discuss the effects of under-sampling and advantages of using parallel imaging that allows faster acquisition [3] .

CS MRI reconstruction problems are extremely computationally intensive and would have taken orders of magnitude longer to perform even a couple decades ago. Sparsi- fying transforms that perform well in CS can be a wide variety from orthogonal to redundant to adaptive dictionaries, each with varying degrees of complexity and applicability. Ning et al. investigate the patched-based trained directional wavelets and extend them into the translation invariant domain to enhance MRI image features [4] . Computationally, orthogonal transforms will operate faster than redundant or adaptive dictionaries which can add considerable time to CS simulations and reconstructions. Sadeghi et al. compare redundant dictionary learning algorithms for computational burden and quality [5] . Now many of these experiments and simulations are able to be implemented and performed fairly quickly due to advances in computer technology and resources. Practical experimentation and application as proof of theory underlay engineering progress. The need for comparing other researcher’s software code and algorithms is essential for prototyping ideas. There is a huge body of research in CS for MRI. How does one evaluate what is the best practical technique for a CS MRI product? Advanced simulation software is necessary for trying techniques and making comparisons.

In previous works, the Translation Invariant Wavelet Transform (TIWT) has been shown to perform exceedingly well in CS against other transforms [6] [7] . Various image types were compared and theoretical analysis of sparsifying by frames was demonstrated [8] . Based on the generalized Dictionary Restricted Isometry Property (D-RIP) for CS, the use of redundant sparsifying transforms like the TIWT are allowed. Candès et al. generalize the RIP for redundant dictionaries in CS [9] . Coifman et al. present the value of the TIWT for de-noising images to suppress visual artifacts that might occur when using a typical wavelet [10] . This transform was found to be advantageous for CS because it also eliminates repetitive line pattern artifacts which may occur when applying a standard decimated orthogonal wavelet transform.

Contribution

The focus here contributes a thorough analysis of the TIWT using under-sampling patterns for use in CS with MRI, not previously provided. The Total Variation (TV) minimization vs. the TIWT sparsifying transform are compared in CS reconstruction for IQ improvements. TV is a calculation of the pixel by pixel finite differences in both horizontal and vertical directions across an image. The TIWT is a redundant circularly shifted wavelet calculation that achieves translational invariance, which is lacking in the standard orthogonal wavelet. Various trajectories that cover a wide range of options and possibilities are used as a testing strategy for proving the TIWT is robust for various noise and under-sampling artifacts. The six different under-sampling patterns analyzed are: uniform random, 1D variable density with large and small center, 2D variable density with and without fully sampled center, uniform radial, and logarithmic spiral. The process of developing and comparing these different measurement matrices also involves processing a number of permutations. Through checking several types of applicable under-sampling patterns in simulation, the goal of proving the robustness of this technique is established. This technique can enhance the reconstruction quality of the CS MRI result and therefore enable the use of less data and acquisition time.

Novel software, FastTestCS, is introduced here for the purpose of reducing computation time in CS reconstruction simulation and testing. It is a tool that can be customized for CS reconstructions using different images, sampling patterns, sparse transforms and optimization techniques. It also is written in the compiled language of C++ that can allow for quick runtime of simulations. Additionally, compute performance and time can be measured more in line with an expected fast product implementation. The FastTestCS executable is designed with a full set of versatile arguments to customize all the needed functionality to run several different simulations.

In summary, a comprehensive analysis of the robustness of the TIWT in CS applications is presented through simulations involving six different under-sampling patterns for MRI. Comparisons in reconstruction compute time are made between the Wavelab toolbox in Matlab and Gnu Scientific Library (GSL) in C++ that show a significant time savings factor of 60×. This transform is computationally feasible as it is shown to be only one order of magnitude more computation than the orthogonal wavelet, yet produces higher fidelity image reconstructions. FastTestCS software is introduced and demonstrated as useful for researchers interested in CS simulations using different images, sparsifying transforms, objective functions and under-sampling patterns in a quick development and test environment that can scale to large studies.

The format consists of five main parts. Section 1 and 2 are introductory and background of theory and pertinent literature on under-sampling in CS for MRI. Section 3 details the methods used for the simulations and development of FastTestCS. This includes under-sampling pattern formula and the CS simulation setup. Section 4 presents results from the simulations of TIWT vs TV CS including images and tables of measurements. Finally, Section 5 discusses these results with the use of FastTestCS.

2. Background: Under-Sampling, Sparsity and Optimization in CS for MRI

The three main ingredients of CS are correct under-sampling, sparsity and optimization. Candès et al. and Lustig et al. state a CS requirement of incoherence between the under-sampling domain and the sparse representation domain [1] [2] . This means that under-sampling in k-space must have artifacts that are incoherent in the linear image reconstruction. A primary aspect of this requirement is determining the best under- sampling technique. Using a uniform random under-sampling matrix could seem appropriate, however, k-space packs a majority of discernible frequency content at its origin which trails off at its periphery. Additionally, just sampling individual points of k-space randomly is non-standard for MRI and may not be as efficient. Sampling is costly in time, therefore, determining the shortest sample trajectory path is desired.

For quick adoption of CS in MRI, patterns that follow existing sequences and trajectories with some under-sampling should be used. Both Cartesian and non-Cartesian under-sampling patterns, such as radial and spiral trajectories, have already been investigated for CS in MRI. Lustig et al. utilize a 1D variable density under-sampling [2] . Bigot et al. investigate a random block of sensing vectors [11] . Castets et al. develop an interleaved golden angle spiral for imaging at high magnetic fields [12] . Chan et al. utilize CS with different radial sampling schemes [13] . Polak et al. consider how the performance of CS is bounded when grouped sampling is necessary [14] . For MRI, sampling speed is somewhat constrained to human physiology and pulse sequence design.

Radial and spiral sampling techniques are used for their advantages. For example, radial acquisitions reduce structural and motion artifacts over standard Cartesian acquisition, while keeping the requirement of artifact incoherence. King improves non- Cartesian spiral scanning with an anisotropic field of view where acquistion is fast and reconstruction improves IQ utilizing gridding [15] . Liu et al. develop a circular Cartesian under-sampling strategy for CS that does not require griding [16] . Ersoz et al. also analyze a Cartesian based radial sampling trajectory for MRI [17] . In many of these cases, due to computational simplicity, the Fast Fourier Transform (FFT) is desired, which operates on a Cartesian grid. Griding and inverse griding with interpolation are often needed to transform from a non-Cartesian sample space to a Cartesian space and back again, adding significant computation when repeated for iterative optimization in CS. However, there are spiral and radial linogram Cartesian based trajectories which do not require griding and interpolation which easily make use of the FFT directly. Therefore, careful selection of sampling patterns does have an impact on reconstruction time.

Sparsity is the amount of an image or signal that can be represented with information. When a signal is n-sparse, the signal has p elements, but only has n elements of valuable information, in some cases,. MRI data and images are not typically sparse, however, a sparsifying transform is used which will compress the image making it sparse in the transform coefficient domain and thereby allowing CS.

Due to under-sampling, the signal recovery in Equation (1) is under-determined. The signal is recovered from corrupted or incomplete measurements and a sampling matrix. An error z is to accommodate real world systems with noise. The error must be sufficiently small with respect to the measurement magnitude, where is at most proportional to the noise level [18] .

(1)

The solution is non-linear and exponentially complex, therefore an alternative constraint is utilized. The problem is now convex and can be minimized with quadratic programming, see Equation (2).

(2)

W represents the sparsifying transform and x is the signal approximation. In order to practically implement this sparse recovery problem, an unconstrained formulation is utilized as shown in Equation (3).

(3)

is a regularization parameter for balancing data consistency and sparsity. This can be solved in many ways; in this work, a Non-Linear Conjugate Gradient (NLCG) method is used, as in [2] [19] .

3. Methods

3.1. MRI Under-Sampling Trajectory Simulation

To prove the effective use of the TIWT as a robust technique for CS MRI, several under-sampling patterns are developed and compared with simulation. The goal of an under-sampling pattern for CS is to provide incoherent and minimal noise like artifacts. Developing under-sampling patterns for MRI trajectories that preserve the data at the center and less at the periphery are of great interest, because the sampling domain is k-space and the majority of the signal power is at the center. The design of these trajectories are described next.

The radial under-sampling technique is calculated in a circular pattern, see Figure 1(a) and Equations (4) and (5). Each line has a unique angle a of. Where L is the number of lines, N is the length of one dimension of the sampling space and n represents a point on that line.

(a) (b) (c)(d) (e) (f)

Figure 1. Under-sampling patterns at 25%: (a) Radial with 70 lines; (b) Spiral; (c) Random; (d) 2D variable density; (e) 2D variable density with 25 point fully sampled center radius; (f) 1D variable density with 32 line fully sampled center.

(4)

(5)

The spiral technique presented here is logarithmic, meaning it will grow from the origin as the distance increases, see Figure 1(b) and Equations (6) and (7). Where a and b are constants that control the radius and expansion of the spiral and t is a point along the spiral

(6)

(7)

Several random distributions are utilized, a basic uniform random under-sampling is used as in Figure 1(c). The variable density random under-sampling in Figure 1(d) uses a 2D radial distance x calculation between the origin and any particular point. This distance is then used in a calculation of a random variable based on a 1D Probability Distribution Function (PDF), that takes a standard deviation () and a distance (x), see Equation (8).

(8)

The resulting random distribution matrix can be thresholded to any desired reduction factor %. Additionally, a fully sampled center can be applied to this technique by setting the probability of the points within the center diameter to 100%, see Figure 1(e).

A 1D variable density under-sampling technique with a fully sampled center is seen in Figure 1(f). This is achieved in much the same way as the 2D variable density calculation, using a 1D distance from the origin of k-space and sampling the entire line.

An analysis is performed for the question of whether to reconstruct CS at a low vs high resolution for the case of a Time of Flight (TOF) Maximum Intensity Projection (MIP). The input high resolution MRI TOF k-space dataset is 256 × 224. It is centered by the highest peak in a 512 × 512 matrix, see Figure 3(a). Notice that this dataset is a partial Fourier encoding, where the center of k-space is fully sampled, and a portion of k-space is not sampled. For simulation of CS, this dataset presents a challenge for reconstruction at low resolution, see Figure 3(b). To use this k-space dataset in simu- lation, it would need to be cropped at 256 × 256 having the peak centered. This does impact IQ of high resolution details, due to the fact that some fine details are deleted in the periphery of k-space. For the simulations performed here, shifting of the k-space peak avoids any cropping of data, and thus minimizes this negative impact to IQ.

(a) (b) (c)(d) (e) (f)

Figure 2. CS reconstructions of T1 MRI with 25% samples. (a) Radial; (b) Spiral; (c) Random; (d) 2D variable density; (e) 2D variable density with 25 point fully sampled center radius; (f) 1D variable density with 32 line fully sampled center.

(a) (b) (c) (d)

Figure 3. TOF MRI k-space data (1 receiver) and MIP reconstruction zoomed in (a) High resolution 512 × 512 data; (b) Cropped low resolution 256 × 256; (c) High resolution reconstruction; (d) Cropped low resolution reconstruction.

Additionally, a shifting of the under-sampling pattern to correlate with the center of k-space is also performed. For CS MRI acquisition, the under-sampling needs to be centered at the k-space peak to be effective and not lose important information at the k-space center and periphery.

3.2. “FastTestCS”, A Fast Software Framework for Testing under-Sampling and Sparsifying Transforms in CS

A new CS testing framework called FastTestCS is proposed here. It is fully written in C++, is modular, with simple methods that can be used similar to functions used in Matlab [20] or other object oriented programs. FastTestCS is built upon the under- standing that CS takes considerably more computation and memory than standard cases of MRI reconstruction. FastTestCS addresses the need to have a more efficient compiled programming environment over prototype Matlab software.

The speed of the TIWT was tested by comparing the transform reconstruction wall time difference between the Wavelab Matlab implementation vs. the transform written in C++. Two resolutions of an image were chosen, 256 × 256 and 512 × 512. The test code was developed by isolating the transform setup, data memory transfers, compu- tation and tear down of the method. These timing tests are performed without the overhead of other CS processing. This time measurement includes calculating 10 for- ward and inverse transforms on both real and imaginary parts of the image for a total of 40 calls to the TIWT method.

This number of transform calculations is comparable to the minimum that must be performed in a CS program. However, often the number of transforms are much greater, such is the case for the NLCG with a quadratic line search and Wolfe con- ditions. Many extra iterations must be performed to determine appropriate search dire- ctions and step sizes.

Four different implementations are tested for the standard orthogonal wavelet and the TIWT. The translational invariance is achieved by performing circular shifts of the data. The orthogonal wavelet provided by Matlab as “WAVEREC2” and “WAVEDEC2”, is compared with the GSL C++ implementation “gsl_wavelet2d_nstransform_matrix”. The TIWT, provided by Wavelab as “FWT2_TI” and “IWT2_TI”, is compared with the FastTestCS implementation “TIDWT” and “ITIDWT”.

The FastTestCS approach implements CS and several sparsifying transforms that run faster in C++ than in Matlab. Algorithms can be prototyped for accuracy first in Matlab, and then re-implemented in C++. The simulations are run in C++ in less time and in parallel across multiple cores and machines. Additionally, FastTestCS has image com- parison tools such as Mean Square Error (MSE), as well as other useful image operations like Normalize, ReadImg or Write Data. It is a Microsoft Visual Studio project and can be compiled for Windows or Linux. It allows sharing code for easy algorithm prototyping, comparison and dataset manipulation. It is a portable and customizable simulation tool to benefit CS research.

There are many required permutations of several parameters and reconstruction tests to run in order to prove robustness. It is apparent that efficiency can be improved by parallelizing computations and simulations, compared with running on just one machine with one thread. Since some CS routines could run over a week to a month, it would be impractical to try many different options serially.

In response, FastTestCS executable is designed with very versatile arguments to implement all the needed functionality and parameters to run several different simulations from the command line. Multiple machines are used to easily distribute many simulations. Simple scripts are used to tar up the software and unique calling parameters, send the package to the specified machine, build and run the package, save the results and transfer the images and results back to a single archiving repository for review. With this method, there are no limits to the number of machines that can be used and simulations that can be achieved more quickly in parallel.

Distribution is achieved for these tests with a small cluster of five individual quad core machines, with one being the source and archive, in addition to having its own portion of processing. Meta-data in the image identifies what simulation parameters produced the image for distinction. Once all the resulting images are together, quantitative and qualitative image comparisons are performed with a “golden” fully sampled image. An executable loads the reconstructed images and the golden image, and makes several calculations for each image such as MSE. The program can then output results to a comma separated table with one line for each image that is easily viewed in a spreadsheet.

MRI k-space data is generated by a General Electric (GE) Healthcare MRI scanner and saved as a “P”-file. GE provided Matlab software to read the P-file and make a k-space matrix. Due to this data being on a Cartesian grid, a fully sampled recon- struction is simply an Inverse FFT. This fully sampled dataset is input to the under- sampling process used in FastTestCS reconstruction by way of a Comma Separated Value (CSV) file. The resulting CS reconstruction is compared to the original fully sampled dataset and measured by MSE to verify its fidelity.

The MSE is defined as a measure of pixel intensity error between and, which are the original and reconstructed image. The smaller this value is, roughly correlates to less measured error and better IQ. MSE is used extensively in this work as a consistent numeric measure of image differences. However, it is important to note that it does not directly correlate to a clinical setting of perceived IQ, where a radiologist’s expertise is needed.

(9)

Three image types are used in the analysis. A phantom image, see Figure 4(a), is a Shepp-Logan phantom, created by Matlab where a 2D FFT converts the image to k-space and saves the CSV. A T1 MRI image, see Figure 4(b), is a single channel T1 MRI brain image where k-space is directly saved to CSV. A TOF MIP Image is a multi slice 8 channel 3D volume TOF brain MRI k-space dataset where each slice and channel image is reconstructed individually and used to generate the MIP, as seen in Figure 4(c) and Figure 4(d). Multiple channels in images are combined using sum of squares.

An important way to analyze different parameters and IQ results of CS recon- structions is to keep all parameters the same in the CS algorithm and just change the transform or sampling pattern. The analysis and results here are done with individual transforms without a combination of other transforms. This is to measure the error of each sampling pattern independently and identify the transform quality.

The steps of the CS reconstruction simulation are as follows and in Figure 5:

1) Start with fully sampled k-space data.

2) Sample the k-space data using a desired under-sampling pattern.

3) Perform the iterative non-linear optimization technique.

4) Compare the error between the fully sampled image and the CS reconstruction to measure the quality of the reconstruction.

The CS reconstruction optimization in FastTestCS is the NLCG with backtracking

(a) (b) (c) (d)

Figure 4. Fully sampled reconstruction simulation images: (a) Shepp-Logan phantom 256 × 256; (b) T1 MRI image 256 × 256; (c) TOF MRI MIP image 192 × 512 × 512 8 channel; (d) TOF MRI (one slice) image 512 × 512.

Figure 5. Implementation of CS image reconstruction simulation.

line search and Wolfe conditions. The formulation of the objective function in Equation (3) requires the calculation of the gradient (), which Lustig, et al. [2] define as Equation (10).

(10)

where can be approximated with the sign function as in Equation (11). The star (*) symbol represents the complex conjugate transpose.

(11)

This unconstrained formulation is versatile when dealing with non-orthonormal transforms. The calculation of the gradients in the unconstrained approach requires the inverse of the transform. For an orthonormal transform, that is simply a multiplication of the coefficients by a transpose of the wavelet dictionary. However, for a non- orthonormal transform, like the TIWT and many of the other transforms, a separate calculation of the inverse transform is used rather than a transpose of the dictionary. In FastTestCS, a transform method is created that calls the forward or inverse of the transform. All coding for the transform is easily integrated, segregated and tested. Real and imaginary parts of MRI data are treated separately throughout all of the optimization.

4. Results

Performance measurements show a big speed up in computation time with optimized C++ code vs Wavelab. See Table 1, where the orthogonal wavelet is optimized in a Matlab toolbox vs in C++ using GSL [22] in FastTestCS, the computational time is very similar. However, when using the TIWT from Wavelab in Matlab vs C++ where the TIWT implementation is coded for FastTestCS, the difference is quite dramatic. FastTestCS is 60 times faster! This shows that coding some computationally time consuming algorithms in C++ would be most beneficial.

Tables 2-8, show a comparison of TIWT vs TV CS with several sampling patterns for IQ analysis on three sample images from Figure 4. These images are diverse in several practical qualities. The Shepp-Logan phantom is commonly used in analysis and

Table 1. Time measurement (seconds) comparison of 40 wavelet transforms (10 complex, forward and reverse), O = Orthogonal, with C++ and Matlab implementations.

Table 2. Radial under-sampling pattern IQ comparison with MSE between CS reconstructed images and reduction factors (%) with 10, 100 and 500 NLCG iterations.

Table 3. Spiral under-sampling pattern IQ comparison with MSE between images and reduction factors.

demonstrates low noise and simple structures. The T1 and TOF MRI brain images portray practical medical images which have differences in noise and fine structure throughout. The 3D TOF additionally represents a sparse image domain which should

Table 4. Random under-sampling pattern IQ comparison with MSE between images and reduction factors.

Table 5. 1D variable density with small (32 line) center under-sampling pattern IQ comparison with MSE between images and reduction factors.

Table 6. 1D variable density with large (60 line) center under-sampling pattern IQ comparison with MSE between images and reduction factors.

Table 7. 2D variable density under-sampling pattern IQ comparison with MSE between images and reduction factors.

lend well to CS due to the use of TV and wavelet minimization functions. Figure 2 portrays the effect each under-sampling pattern has on the linear reconstruction of Figure 4(b).

Table 8. 2D variable density with small (25 pixel radius) center under-sampling pattern IQ comparison with MSE between images and reduction factors.

In Table 2, there is a linear vs CS reconstruction IQ comparison with a variable density radial under-sampling pattern across the three images. There is also an analysis of stopping criteria where MSE results of 10 iterations are compared with 100 and 500 iterations. Note for the phantom, the CS reconstruction with the TIWT has a much lower MSE than that of TV at 10 NLCG iterations and at a lower under-sampling rate of 25%. However, at 100 and 500 iterations the MSE is similar. This shows that 10 iterations is insufficient and performing more than 100 iterations may not change the image substantially. Through this data, there also appears to be different iteration limits to improving IQ for certain image types and transforms. When looking for a robust general solution, this experiment highlights how important it is to have an analysis to determine good iteration stopping criteria. Given this challenge, in order to have a consistent analysis, 100 iterations are chosen as the stopping criteria for all the other simulations.

In Table 3 for spiral under-sampling, comparisons can be made with radial sampling in terms of MSE. For instance, in the TOF MIP at 38% samples for TIWT and TV, the MSE is 10.4 and 33.0, respectively, which is similar to the results shown in Table 2. Looking at the same under-sampling ratios, a better IQ is achieved in an initial linear reconstruction. However, this does not always correlate to a better CS reconstruction, as seen in some results of this table.

Table 4 shows results using a uniform random under-sampling pattern. Much higher MSE is observed in all linear and CS reconstructions. These results show it is important to have dense sampling at the center of k-space for reconstruction. Notice the range of MSE at lower % vs high %, even at 75% samples, the MSE of the recon- struction is similar with radial under-sampling at 25%. However, even with this poor choice of under-sampling, the TIWT does operate effectively.

In Table 4, according to the MSE of the TOF MIP, it appears that the IQ does not improve very much when CS is utilized, and in some cases for TV, even degrades. This highlights the problem using a uniform random under-sampling in MRI. The majority of large structure information is at the center of k-space, and some of that information is not being recovered which can drastically change reconstruction. Although, this analysis is not entirely negative, as seen by inspecting intermediate CS reconstructed TOF images that make up the MIP. Figure 6(b) and Figure 7(a) show a 38% under- sampling linear reconstruction of the TOF slice and MIP. A fully sampled slice is shown in Figure 6(a). Intense noise has clouded the TOF and MIP image. After a TIWT or TV

(a) (b) (c) (d)

Figure 6. TOF reconstruction with 38% random under-sampling, receiver 1, slice 1. (a) Fully sampled; (b) Linear; (c) TIWT CS; (d) TV CS.

(a) (b) (c)

Figure 7. TOF MIP reconstruction with 38% random under-sampling. (a) Linear; (b) TIWT CS; (c) TV CS.

CS reconstruction, the noise is significantly reduced, see Figure 6(c), Figure 7(b), Figure 6(d) and Figure 7(c). Notice that comparing these results and the fully sampled MIP, Figure 4(c), shows a large difference in background intensity. This is greatly reduced in the CS reconstruction and may be a desired trait because it intensifies the Signal to Noise Ratio (SNR). It is important to note that some fine details and other large structure intensities are reconstructed slightly different between TV and TIWT. In terms of the robustness in MSE of the TIWT, it outperforms TV in all cases.

Table 5 and Table 6 represent the use of several 1D variable density under- samplings of k-space for good IQ and simple adoption into standard Cartesian MRI pulse sequences. In Table 6, a larger center of 2D fully sampled k-space contributes to a much better linear reconstruction, which in turn, produces an improved CS recon- struction over the smaller center shown in Table 5. Additionally, these consistent improvements do bolster a robust use of TIWT for CS.

Table 7 shows the use of a random variable density 2D under-sampling. Similar to the fully uniform random under-sampling cases of Table 4, there are drastic differences in MSE for reconstructions at different under-sampling percentages. Lower sampling results in less improvement in the CS reconstructions, however, this result improves when more samples are taken. Notice that for both these cases when a low percentage of samples are taken, the center of k-space is not fully sampled and the image recon- structions suffer a large loss of quality in major image structures and intensity. It appears that this under-sampling pattern should also be avoided.

Table 8 shows improved results over Table 7 by modifying the 2D under-sampling pattern with the addition of a center that is fully sampled. These results are very promising and show that this under-sampling pattern produces some of the best reconstructions. The TIWT also performs quite well with this sampling pattern.

Figure 3(c) and Figure 3(d) represent an enlarged portion of the fully sampled TOF MIP. A high resolution (512 × 512) linear reconstruction is compared with a cropped k-space low resolution (256 × 256) reconstruction zipped to 512 × 512 for MIP calculation, The MSE difference between the two results is about 3.0, which is small compared to some under-sampled linear reconstructions. However, in a visual analysis of the images, some very fine vessels are missing or blurred in the low resolution reconstruction causing an IQ impairment.

To mitigate this, the MIP simulations shown in the tables are all done by shifting the k-space data so that it is not cropped, but rather has the center of k-space and the sampling center match slightly off center. The MSE difference in this case is 0.11, which is negligible, and verifies the assumption that there is no visible difference in IQ between low and high resolution fully sampled reconstructions. With this accuracy, there is confidence in running these tests at low resolution with a shifted k-space.

Performing all reconstructions at a high resolution would come at a high compu- tation cost of about 5× longer, therefore, lower resolution reconstructions are perfor- med. As a time comparison, a 512 × 512 8 channel 192 slice MIP CS reconstruction with the TIWT for 10 iterations takes 20 hours, whereas, the same reconstruction at 256 × 256 takes 4 hours, single threaded, with an Intel Core i7 computer.

5. Discussion

5.1. Good under-Sampling Patterns

Here it is shown that the choice of good under-sampling patterns clearly has a correlation to the data being gathered. The center of k-space visually has a more intense magnitude of data, therefore, sampling density at the center contributes to a better SNR and a lower MSE. The radial, spiral, 1D and 2D variable density with the center fully sampled performs good, compared to the other patterns that did not sample with a dense center. The balancing of image quality with how the under-sampling pattern dissipates, is a complex analysis that requires visual as well as error analysis.

5.2. CS Reconstruction for MRI

As seen in the synthetic phantom error tables, some CS reconstructions may perform exceptionately well in comparison to real in-vivo MRI data. The lack of noise and the few distinct edges in the original phantom data contribute to a high success for the TIWT and especially TV, which measures differences between pixels. With that in mind, care must be taken for each MRI imaging application that uses CS due to likely differences in noise statistics and image types affecting the IQ performance of the CS reconstruction. This hightlights the value of having a broad spectrum of tests so that the correct CS parameters can be chosen and used.

5.3. Reconstruction of MIP

CS reconstruction of TOF data used in a MIP can be performed at a lower acquisition resolution and then transformed to a higher resolution because there are no visible differences in smaller vessels. The analysis of IQ with MSE is challenging for MIP because the error measurement may not pick up the changes in small vessels being of little numerical error. Additionally, changes in the background intensity of the image do contribute to larger error. Having a larger MSE may be misleading compared to what is observed with a visual analysis.

Special care must be taken for under-sampling with partially encoded Fourier k-space data. Visually, the image generated without small amounts of data at the periphery does have a negative impact on small vessel details of the MIP CS reconstruction. Therefore, adjusting the k-space center is needed to prevent cropping any data and have more meaningful IQ comparisons.

5.4. Stopping Criteria

The stopping criteria is determined by analyzing the IQ after a certain number of iterations. This proves to be a very challenging issue due to the fact that different image types will require a different number of iterations to produce the best reconstructions. Image types with very little noise require many more iterations, however, the recon- struction IQ improves greatly as well. This improvement is probably due to the fundamental aspect of CS reconstruction theory of sparsity. When higher sparsity can be achieved in transforming these images, CS has a greater probability of better recon- structions with less data. However, most real images do have noise.

Despite this dilemma, a solution can be observed where enough iterations can be achieved that produces a result that will not change IQ significantly if more iterations are performed. Then this number can be used for other similar tests. This aspect of CS highlights the importance of performing an analysis for determining the best stopping criteria for specific image types and algorithms.

5.5. Computation Time with FastTestCS

The choice to develop FastTestCS came out of the necessity to speed up reconstructions using a more complex TIWT in CS with larger 3D data sets. Table 1 proves a conside- rable speedup is achieved by coding in a compiled language like C++ which is a more realistic time measurement of reconstruction than using less optimized software. The use of the TIWT vs a standard orthogonal wavelet is 10× more processing, achieving an improvment of reconstruction IQ. Although, this is at a high compute burden, it is less than more complex algorithms that incorporate dictionary learning. Having faster software contributes to more research opportunities such as being able to compare more sampling patterns, finding better stopping criteria, and performing higher resolution experiments with more images. For instance, comparing the CS recon- struction time of the TOF MIP, which uses +1500 images (256 × 256), the Wavelab TIWT implementation in Matlab would have taken approximately 10 days to recon- struct, whereas, with FastTestCS, this same reconstruction takes 4 hours on a single compute core. Additionally, by adding parallelization, with five quad-core machines, this reconstruction takes less than 15 minutes. There is also great value in FastTestCS to assess true performance of sparsifying transforms because it is better to compare a fast compiled language implementation vs a non-optimized version. This is the case in the use of the TIWT from Wavelab resulting in a 60x longer compute time compared with FastTextCS.

5.6. Use of FastTestCS

FastTestCS may be customized to include any image or signal, objective function, sampling pattern and sparsifying transform to test different scenarios and compare reconstructions. Additionally, any libraries or algorithms available in C++ can be incorporated into FastTestCS. Having consistent CS input data and algorithms to test various options provides a clear analysis and simulation tool with unlimited possibilities for future enhancements and tests. The executable and command line interface provides a convenient way to parallelize on different compute cores and machines to quickly test across various inputs and options in an automated way. Contact the author for software availability to obtain reproducible results, usage and advancement, http://people.uwm.edu/gxu4uwm/fasttestcs.

5.7. Future Work

The applicability of using the TIWT for CS MRI is wide spread due to the variety of image types and sampling patterns encountered in MRI. Therefore, there are many applications along those lines that could be researched in the future. Additionaly, different wavelet filters and settings could be investigated to find which applications perform best with which filters to produce better IQ. Also, future expansions to FastTestCS are possible to add functionaltiy for testing other images, sparse transforms, sampling patterns and objective functions.

6. Conclusions

By simulation of CS reconstructions and measurement of IQ differences, results indicate robustness of the TIWT. Several comprehensive under-sampling patterns are designed to compare the TIWT with TV as sparsifying transforms in CS MRI recon- struction. The TIWT performs consistently well in terms of MSE for all under- sampling patterns. An in-depth look at high and low resolution TOF MIP CS recon- struction is evaluated for IQ differences. Analysis shows that careful under-sampling must be performed at the center of k-space to preserve important data for better recon- struction. With the TIWT, good k-space under-sampling and a reliable CS recon- struction algorithm, improved IQ and greater speed-up for CS MRI can be achieved.

FastTestCS software is introduced and demonstrated as a helpful tool for CS reconstruction research that offers the needed speedup with process distribution to perform experiments with large datasets. In the case of the TIWT, a 60× time improve- ment was achieved. FastTestCS also is easily customized with each under-sampling pattern used in this investigation.

Acknowledgements

The author wishes to thank Professor Guangwu Xu for guidance, Academic Advisor, Department of Computer Science and Engineering at the University of Wisconsin- Milwaukee, and Dr. Kevin King, Scientist at GE Healthcare for providing the MRI data along with helpful review.

References

[1] Candès, Emmanuel, J. and Wakin, M.B. (2008) An Introduction to Compressive Sampling. IEEE Signal Processing Magazine, 25, 21-30.

http://dx.doi.org/10.1109/MSP.2007.914731

[2] Lustig, M., Donoho, D. and Pauly, J.M. (2007) Sparse MRI: The Application of Compressed Sensing for Rapid MR Imaging. Magnetic Resonance in Medicine, 58, 1182-1195.

http://dx.doi.org/10.1002/mrm.21391

[3] Deshmane, A., Gulani, V., Griswold, M.A. and Seiberlich, N. (2012) Parallel MR Imaging. Journal of Magnetic Resonance Imaging, 36, 55-72.

http://dx.doi.org/10.1002/jmri.23639

[4] Ning, B., Qu, X., Guo, D., Hu, C. and Chen, Z. (2013) Magnetic Resonance Image Reconstruction Using Trained Geometric Directions in 2D Redundant Wavelets Domain and Non-Convex Optimization. Magnetic Resonance Imaging, 31, 1611-1622.

http://dx.doi.org/10.1016/j.mri.2013.07.010

[5] Sadeghi, M., Babaie-Zadeh, M. and Jutten, C. (2013) Learning Overcomplete Dictionaries Based on Parallel Atom-Updating. 2013 IEEE International Workshop on Machine Learning for Signal Processing (MLSP), Southampton, 22-25 September 2013, 1-5.

http://dx.doi.org/10.1109/MLSP.2013.6661975

[6] Baker, C.A., King, K.F., Liang, D. and Ying, L. (2009) Investigation of Sparsifying Transforms for Compressed Sensing in MR Reconstruction. 17th ISMRM Annual Meeting, Honolulu, 18-24 April 2009, 4583.

[7] Baker, C.A., King, K.F., Liang, D. and Ying, L. (2011) Translational-Invariant Dictionaries for Compressed Sensing in Magnetic Resonance Imaging. 2011 IEEE International Symposium on Biomedical Imaging: From Nano to Macro, Chicago, 30 March-2 April 2011, 1602-1605.

http://dx.doi.org/10.1109/ISBI.2011.5872709

[8] Baker, C.A. (2016) Sparse Representation by Frames with Signal Analysis. Journal of Signal and Information Processing, 7, 39-48.

http://dx.doi.org/10.4236/jsip.2016.71006

[9] Candès, E.J., Eldar, Y.C., Needell, D. and Randall, P. (2011) Compressed Sensing with Coherent and Redundant Dictionaries. Applied and Computational Harmonic Analysis, 31, 59-73.

http://dx.doi.org/10.1016/j.acha.2010.10.002

[10] Coifman, R.R. and Donoho, D.L. (1995) Translation-Invariant De-Noising. Springer, New York.

http://dx.doi.org/10.1007/978-1-4612-2544-7_9

[11] Bigot, J., Boyer, C. and Weiss, P. (2016) An Analysis of Block Sampling Strategies in Compressed Sensing. IEEE Transactions on Information Theory, 62, 2125-2139.

http://dx.doi.org/10.1109/TIT.2016.2524628

[12] Castets, C.R., Lefrancois, W., Wecker, D., Ribot, E.J., Trotier, A.J., Thiaudière, E., Franconi, J.M. and Miraux, S. (2016) Fast 3D Ultrashort Echo-Time Spiral Projection Imaging Using Golden-Angle: A Flexible Protocol for in Vivo Mouse Imaging at High Magnetic Field. Magnetic Resonance in Medicine, Epub.

[13] Chan, R.W., Ramsay, E.A., Cheung, E.Y. and Plewes, D.B. (2012) The Influence of Radial Undersampling Schemes on Compressed Sensing Reconstruction in Breast MRI. Magnetic Resonance in Medicine, 67, 363-377.

http://dx.doi.org/10.1002/mrm.23008

[14] Polak, A.C., Duarte, M.F. and Goeckel, D.L. (2015) Performance Bounds for Grouped Incoherent Measurements in Compressive Sensing. IEEE Transactions on Signal Processing, 63, 2877-2887.

http://dx.doi.org/10.1109/TSP.2015.2412912

[15] King, K.F. (1998) Spiral Scanning with Anisotropic Field of view. Magnetic Resonance in Medicine, 39, 448-456.

http://dx.doi.org/10.1002/mrm.1910390315

[16] Liu, J. and Saloner, D. (2014) Accelerated MRI with CIRcular Cartesian under Sampling (CIRCUS): A Variable Density Cartesian Sampling Strategy for Compressed Sensing and Parallel Imaging. Quantitative Imaging in Medicine and Surgery, 4, 57-67.

[17] Ersoz, A. and Tugan Muftuler, L. (2016) Undersampled Linogram Trajectory for Fast Imaging (ULTI): Experiments at 3 T and 7 T. NMR in Biomedicine, 29, 340-348.

http://dx.doi.org/10.1002/nbm.3475

[18] Candès, E.J., Romberg, J.K. and Tao, T. (2006) Stable Signal Recovery from Incomplete and Inaccurate Measurements. Communications on Pure and Applied Mathematics, 59, 1207-1223.

http://dx.doi.org/10.1002/cpa.20124

[19] Nocedal, J. and Wright, S. (2006) Numerical Optimization. Springer, Berlin.

[20] MathWorks Inc. (2005) MATLAB: The Language of Technical Computing. Desktop Tools and Development Environment. MathWorks, Version 7, Vol. 9.

[21] Donoho, D., Maleki, A. and Shahram, M. (2006) Wavelab 850. Software Toolkit for Time-Frequency Analysis.

[22] Gough, B. (2009) GNU Scientific Library Reference Manual. Network Theory Ltd.