This example demonstrates solving a convolutional sparse coding problem with a greyscale dictionary and a colour signal [36]
$$\mathrm{argmin}_\mathbf{x} \; (1/2) \sum_c \left\| \sum_m \mathbf{d}_m * \mathbf{x}_{c,m} - \mathbf{s}_c \right\|_2^2 + \lambda \sum_c \sum_m \| \mathbf{x}_{c,m} \|_1 \;,$$
where $\mathbf{d}_{m}$ is the $m^{\text{th}}$ dictionary filter, $\mathbf{x}_{c,m}$ is the coefficient map corresponding to the $m^{\text{th}}$ dictionary filter and channel $c$ of the input image, and $\mathbf{s}_c$ is channel $c$ of the input image.
from __future__ import print_function
from builtins import input
from builtins import range
import pyfftw # See https://github.com/pyFFTW/pyFFTW/issues/40
import numpy as np
from sporco import util
from sporco import plot
plot.config_notebook_plotting()
import sporco.metric as sm
from sporco.admm import cbpdn
Load example image.
img = util.ExampleImages().image('kodim23.png', scaled=True,
idxexp=np.s_[160:416,60:316])
Highpass filter example image.
npd = 16
fltlmbd = 10
sl, sh = util.tikhonov_filter(img, fltlmbd, npd)
Load greyscale dictionary and display it.
D = util.convdicts()['G:8x8x64']
plot.imview(util.tiledict(D), fgsz=(7, 7))
Set admm.cbpdn.ConvBPDN solver options.
lmbda = 1e-1
opt = cbpdn.ConvBPDN.Options({'Verbose': True, 'MaxMainIter': 200,
'RelStopTol': 5e-3, 'AuxVarObj': False})
Initialise and run CSC solver.
b = cbpdn.ConvBPDN(D, sh, lmbda, opt, dimK=0)
X = b.solve()
print("ConvBPDN solve time: %.2fs" % b.timer.elapsed('solve'))
Itn Fnc DFid Regℓ1 r s ρ ---------------------------------------------------------------- 0 4.97e+02 3.05e+00 4.94e+03 9.92e-01 3.61e-02 6.00e+00 1 4.24e+02 1.66e+01 4.07e+03 8.93e-01 1.64e-01 6.00e+00 2 3.46e+02 3.16e+01 3.14e+03 4.46e-01 2.84e-01 9.89e+00 3 3.43e+02 3.79e+01 3.05e+03 2.77e-01 2.38e-01 8.76e+00 4 3.42e+02 3.67e+01 3.06e+03 2.38e-01 1.49e-01 6.69e+00 5 3.03e+02 3.63e+01 2.66e+03 2.08e-01 1.17e-01 5.98e+00 6 2.71e+02 3.76e+01 2.33e+03 1.77e-01 8.89e-02 5.98e+00 7 2.49e+02 3.96e+01 2.10e+03 1.46e-01 6.93e-02 5.98e+00 8 2.36e+02 4.15e+01 1.94e+03 1.16e-01 6.00e-02 5.98e+00 9 2.30e+02 4.26e+01 1.87e+03 9.43e-02 5.32e-02 5.98e+00 10 2.25e+02 4.32e+01 1.82e+03 7.92e-02 4.63e-02 5.98e+00 11 2.19e+02 4.33e+01 1.75e+03 6.82e-02 4.14e-02 5.98e+00 12 2.11e+02 4.35e+01 1.68e+03 6.28e-02 3.78e-02 5.43e+00 13 2.06e+02 4.37e+01 1.62e+03 5.88e-02 3.31e-02 4.95e+00 14 2.02e+02 4.41e+01 1.58e+03 5.24e-02 2.88e-02 4.95e+00 15 1.97e+02 4.44e+01 1.53e+03 4.64e-02 2.61e-02 4.95e+00 16 1.93e+02 4.48e+01 1.49e+03 4.10e-02 2.46e-02 4.95e+00 17 1.91e+02 4.50e+01 1.46e+03 3.68e-02 2.29e-02 4.95e+00 18 1.89e+02 4.51e+01 1.44e+03 3.54e-02 2.09e-02 4.43e+00 19 1.88e+02 4.52e+01 1.43e+03 3.25e-02 1.90e-02 4.43e+00 20 1.86e+02 4.52e+01 1.40e+03 2.98e-02 1.77e-02 4.43e+00 21 1.83e+02 4.53e+01 1.37e+03 2.74e-02 1.67e-02 4.43e+00 22 1.80e+02 4.54e+01 1.35e+03 2.67e-02 1.56e-02 4.02e+00 23 1.79e+02 4.55e+01 1.33e+03 2.49e-02 1.43e-02 4.02e+00 24 1.77e+02 4.57e+01 1.32e+03 2.32e-02 1.32e-02 4.02e+00 25 1.76e+02 4.58e+01 1.30e+03 2.16e-02 1.24e-02 4.02e+00 26 1.75e+02 4.59e+01 1.29e+03 2.02e-02 1.18e-02 4.02e+00 27 1.73e+02 4.59e+01 1.28e+03 1.89e-02 1.12e-02 4.02e+00 28 1.72e+02 4.60e+01 1.26e+03 1.78e-02 1.07e-02 4.02e+00 29 1.71e+02 4.60e+01 1.25e+03 1.75e-02 1.01e-02 3.67e+00 30 1.70e+02 4.61e+01 1.24e+03 1.66e-02 9.60e-03 3.67e+00 31 1.69e+02 4.62e+01 1.23e+03 1.58e-02 9.06e-03 3.67e+00 32 1.68e+02 4.62e+01 1.22e+03 1.50e-02 8.58e-03 3.67e+00 33 1.67e+02 4.63e+01 1.21e+03 1.43e-02 8.17e-03 3.67e+00 34 1.66e+02 4.64e+01 1.20e+03 1.35e-02 7.80e-03 3.67e+00 35 1.66e+02 4.64e+01 1.19e+03 1.29e-02 7.44e-03 3.67e+00 36 1.65e+02 4.65e+01 1.18e+03 1.23e-02 7.10e-03 3.67e+00 37 1.64e+02 4.65e+01 1.18e+03 1.17e-02 6.79e-03 3.67e+00 38 1.64e+02 4.66e+01 1.17e+03 1.12e-02 6.53e-03 3.67e+00 39 1.63e+02 4.66e+01 1.16e+03 1.06e-02 6.30e-03 3.67e+00 40 1.62e+02 4.66e+01 1.16e+03 1.02e-02 6.06e-03 3.67e+00 41 1.62e+02 4.67e+01 1.15e+03 9.73e-03 5.78e-03 3.67e+00 42 1.61e+02 4.67e+01 1.14e+03 9.33e-03 5.53e-03 3.67e+00 43 1.60e+02 4.68e+01 1.14e+03 8.93e-03 5.30e-03 3.67e+00 44 1.60e+02 4.68e+01 1.13e+03 8.55e-03 5.11e-03 3.67e+00 45 1.59e+02 4.68e+01 1.12e+03 8.20e-03 4.94e-03 3.67e+00 46 1.59e+02 4.69e+01 1.12e+03 8.24e-03 4.76e-03 3.34e+00 47 1.59e+02 4.69e+01 1.12e+03 7.96e-03 4.55e-03 3.34e+00 48 1.58e+02 4.69e+01 1.11e+03 7.68e-03 4.35e-03 3.34e+00 49 1.58e+02 4.70e+01 1.11e+03 7.40e-03 4.18e-03 3.34e+00 50 1.58e+02 4.70e+01 1.11e+03 7.12e-03 4.03e-03 3.34e+00 51 1.57e+02 4.70e+01 1.10e+03 6.85e-03 3.89e-03 3.34e+00 52 1.57e+02 4.70e+01 1.10e+03 6.60e-03 3.75e-03 3.34e+00 53 1.57e+02 4.70e+01 1.10e+03 6.36e-03 3.61e-03 3.34e+00 54 1.56e+02 4.71e+01 1.09e+03 6.13e-03 3.48e-03 3.34e+00 55 1.56e+02 4.71e+01 1.09e+03 5.91e-03 3.36e-03 3.34e+00 56 1.56e+02 4.71e+01 1.08e+03 5.70e-03 3.25e-03 3.34e+00 57 1.55e+02 4.71e+01 1.08e+03 5.49e-03 3.15e-03 3.34e+00 58 1.55e+02 4.72e+01 1.08e+03 5.30e-03 3.05e-03 3.34e+00 59 1.55e+02 4.72e+01 1.07e+03 5.12e-03 2.96e-03 3.34e+00 60 1.54e+02 4.72e+01 1.07e+03 4.94e-03 2.86e-03 3.34e+00 ---------------------------------------------------------------- ConvBPDN solve time: 66.44s
Reconstruct image from sparse representation.
shr = b.reconstruct().squeeze()
imgr = sl + shr
print("Reconstruction PSNR: %.2fdB\n" % sm.psnr(img, imgr))
Reconstruction PSNR: 33.13dB
Display low pass component and sum of absolute values of coefficient maps of highpass component.
gamma = lambda x, g: x**g
fig = plot.figure(figsize=(14, 7))
plot.subplot(1, 2, 1)
plot.imview(sl, title='Lowpass component', fig=fig)
plot.subplot(1, 2, 2)
plot.imview(gamma(np.sum(abs(X), axis=b.cri.axisM).squeeze(), 0.4),
title='Sparse representation', fig=fig)
fig.show()
Display original and reconstructed images.
fig = plot.figure(figsize=(14, 7))
plot.subplot(1, 2, 1)
plot.imview(img, title='Original', fig=fig)
plot.subplot(1, 2, 2)
plot.imview(imgr, title='Reconstructed', fig=fig)
fig.show()
Get iterations statistics from solver object and plot functional value, ADMM primary and dual residuals, and automatically adjusted ADMM penalty parameter against the iteration number.
its = b.getitstat()
fig = plot.figure(figsize=(20, 5))
plot.subplot(1, 3, 1)
plot.plot(its.ObjFun, xlbl='Iterations', ylbl='Functional', fig=fig)
plot.subplot(1, 3, 2)
plot.plot(np.vstack((its.PrimalRsdl, its.DualRsdl)).T,
ptyp='semilogy', xlbl='Iterations', ylbl='Residual',
lgnd=['Primal', 'Dual'], fig=fig)
plot.subplot(1, 3, 3)
plot.plot(its.Rho, xlbl='Iterations', ylbl='Penalty Parameter', fig=fig)
fig.show()