SPIRAL: Sparse Poisson Intensity Reconstruction ALgorithms [Download]

Overview
The Sparse Poisson Intensity Reconstruction ALgrotihms (SPIRAL) toolbox, SPIRALTAP.m, is MATLAB code for recovering sparse signals from Poisson observations. SPIRAL minimizes a regularized negative log-likelihood objective function with various penalty choices for the regularization terms:
  • Sparsity (l1 norm) of the coefficients in an orthonormal basis [SPIRAL-ONB],
  • Total variation seminorm of the image [SPIRAL-TV],
  • Penalty based on Recursive Dyadic Partitions (RDPs) [SPIRAL-RDP], and
  • Penalty based on translationally-invariant (cycle-spun) RDPs [SPIRAL-RDPTI].
For more details, see: Zachary T. Harmany, Roummel F. Marcia, Rebecca M. Willett, “This is SPIRAL-TAP: Sparse Poisson Intensity Reconstruction ALgorithms – Theory and Practice,” Submitted to IEEE Transactions on Image Processing. [PDF (arXiv.org)]

To aid users we provide a few examples of our algorithm. To view a demonstration, execute in MATLAB
>> SPIRALdemo
the varialbe demo in SPIRALdemo.m to 1, 2, or 3 selects among three simulations. Details of each can be found below.

NOTE: Some of these demonstrations utilize the Rice Wavelet Toolbox (freely available online) to compute the discrete wavelet transform. We include this toolbox here (although it may need to be recompiled on your platform). We also use the FISTA algorithm (in the denoise directory) of Beck and Teboulle for constrained Total Variation denoising.

Demonstration 1
Description: One dimensional compressed sensing example penalizing the sparsity (l1 norm) of the coefficients in the canonical basis. Here the true signal f is of length 100,000 with 1,500 nonzero entries yielding a sparsity of 1.5%. We take 40,000 compressive measurements in y. The average number of photons per measurement is 15.03, with a maximum of 145. We run SPIRAL until the relative change in the iterates falls below a tolerance of 1x10^-8, up to a maximum of 100 iterations (however only 37 iterations are required to satisfy the stopping criterion).

Output: This demonstration automatically displays the following:




Demo1_Setup
Figure 1: Simulation setup (true signal, true detector intensity, observed counts) Note: Zoomed to first 200 samples only.




Demo1_fhatSPIRAL
Figure 2: Reconstructed signal (red) overlaid on top of the true signal (blue) Note: Zoomed to first 200 samples only.




Demo1_RMSEEvolution
Figure 3: RMSE error evolution versus compute time




Demo1_Objective
Figure 4: Objective evolution versus compute time



Demonstration 2
Description: Here we consider an image deblurring example. The true signal is a 128x128 Shepp-Logan phantom image with mean intensity 1.22x10^5. The true detector mean intensity is 45.8, and the observed photon count mean is 45.8 with a maximum of 398. Here we consider four penalization methods:
  • Sparsity (l1 norm) of coefficients in an orthonormal (wavelet) basis,
  • Total variation of the image,
  • Penalty based on Recursive Dyadic Partitions (RDPs), and
  • Penalty based on Translationally-Invariant (cycle-spun) RDPs.
We run all the SPIRAL methods for a minimum of 50 iterations until the relative change in the iterates falls below a tolerance of 1x10^-8, up to a maximum of 100 iterations (however only ~70 iterations are required to satisfy the stopping criterion for all the methods).

Output: This demonstration automatically displays the following [since this is an image processing example, we explicitly show the outputs]:

True SignalTrue Detector IntensityObserved Photon Counts
Demo2_fDemo2_AfDemo2_y
Figure 1: Simulation setup




Demo2_Objective
Figure 2: The objective evolution for the methods where explicit computation of the objective is possible




Demo2_RMSEEvolution
Figure 3: Normalized RMS error evolution versus computation time



SPIRAL-ONBSPIRAL-TV
RMSE (%) = 24.83RMSE (%) = 25.97

Demo2_fhatSPIRALonb

Demo2_fhatSPIRALtv

SPIRAL-RDPSPIRAL-RDP-TI
RMSE (%) = 24.83RMSE (%) = 25.97

Demo2_fhatSPIRALrdp

Demo2_fhatSPIRALrdpti

Figure 4: Final reconstructions



SPIRAL-ONBSPIRAL-TV
RMSE (%) = 24.83RMSE (%) = 25.97

Demo2_diffSPIRALonb Demo2_Colorbar

Demo2_diffSPIRALtv Demo2_Colorbar

SPIRAL-RDPSPIRAL-RDP-TI
RMSE (%) = 24.83RMSE (%) = 25.97

Demo2_diffSPIRALrdp Demo2_Colorbar

Demo2_diffSPIRALrdpti Demo2_Colorbar

Figure 5: Magnitude of the error between the final reconstruction and the true phantom image. Note: Color scale indicates percentage difference over all reconstructions.


Demonstration 3
Description: This demonstration is similar to Demonstration 2, but uses the 256x256 'Cameraman' image instead of the 128x128 Shepp-Logan phantom. The true signal has a mean intensity of 1.19x10^5, the true detector mean intensity is 44.4, and the observed photon count mean is 44.4 with a maximum of 111. Due to the larger problem size, we need to alter the termination criteria. We run all the SPIRAL methods for a minimum of 50 iterations until the relative change in the iterates falls below a tolerance of 1x10^-8, up to a maximum of 300 iterations.

Output: Like Demonstration 2, this demonstration automatically displays the following:

True SignalTrue Detector IntensityObserved Photon Counts

Demo3_f

Demo3_Af

Demo3_y

Figure 1: Simulation setup




Demo3_Objective
Figure 2: The objective evolution for the methods where explicit computation of the objective is possible




Demo3_RMSEEvolution
Figure 3: Normalized RMS error evolution versus computation time



SPIRAL-ONBSPIRAL-TV
RMSE (%) = 10.92RMSE (%) = 9.88

Demo3_fhatSPIRALonb

Demo3_fhatSPIRALtv

SPIRAL-RDPSPIRAL-RDP-TI
RMSE (%) = 11.44RMSE (%) = 9.99

Demo3_fhatSPIRALrdp

Demo3_fhatSPIRALrdpti

Figure 4: Final reconstructions




SPIRAL-ONBSPIRAL-TV
RMSE (%) = 10.92RMSE (%) = 9.88

Demo3_diffSPIRALonb Demo3_Colorbar

Demo3_diffSPIRALtv Demo3_Colorbar

SPIRAL-RDPSPIRAL-RDP-TI
RMSE (%) = 11.44RMSE (%) = 9.99

Demo3_diffSPIRALrdp Demo3_Colorbar

Demo3_diffSPIRALrdpti Demo3_Colorbar

Figure 5: Magnitude of the error between the final reconstruction and the true phantom image. Note: Color scale indicates percentage difference over all reconstructions.