Colour ℓ1-TV Denoising (CuPy Version)¶

This example demonstrates the use of class tvl1.TVL1Denoise for removing salt & pepper noise from a colour image using Total Variation regularization with an ℓ1 data fidelity term (ℓ1-TV denoising). This variant of the example uses the GPU accelerated version of tvl1 within the sporco.cupy subpackage.

from __future__ import print_function
from builtins import input

import numpy as np

from sporco import util
from sporco import signal
from sporco import metric
from sporco import plot
plot.config_notebook_plotting()
from sporco.cupy import (cupy_enabled, np2cp, cp2np, select_device_by_load,
                         gpu_info)
from sporco.cupy.admm import tvl1

Load reference image.

img = util.ExampleImages().image('monarch.png', scaled=True,
                                 idxexp=np.s_[:,160:672])

Construct test image corrupted by 20% salt & pepper noise.

np.random.seed(12345)
imgn = signal.spnoise(img, 0.2)

Set regularization parameter and options for ℓ1-TV denoising solver. The regularization parameter used here has been manually selected for good performance.

lmbda = 8e-1
opt = tvl1.TVL1Denoise.Options({'Verbose': True, 'MaxMainIter': 200,
                                'RelStopTol': 5e-3, 'gEvalY': False,
                                'AutoRho': {'Enabled': True}})

Create solver object and solve, returning the the denoised image imgr.

if not cupy_enabled():
    print('CuPy/GPU device not available: running without GPU acceleration\n')
else:
    id = select_device_by_load()
    info = gpu_info()
    if info:
        print('Running on GPU %d (%s)\n' % (id, info[id].name))

b = tvl1.TVL1Denoise(np2cp(imgn), lmbda, opt)
imgr = cp2np(b.solve())

Running on GPU 0 (NVIDIA GeForce RTX 2080 Ti)

Itn   Fnc       DFid      RegTV     r         s         ρ
----------------------------------------------------------------
1.67e+05  8.07e+04  1.08e+05  3.60e-01  1.29e+00  1.70e+00
1.43e+05  9.34e+04  6.25e+04  2.65e-01  8.24e-01  1.70e+00
1.92e+05  1.18e+05  9.20e+04  2.88e-01  4.24e-01  9.63e-01
1.89e+05  1.29e+05  7.46e+04  2.67e-01  2.65e-01  7.94e-01
1.59e+05  1.14e+05  5.63e+04  1.92e-01  2.94e-01  7.94e-01
1.41e+05  9.96e+04  5.12e+04  1.47e-01  2.03e-01  6.42e-01
1.30e+05  9.28e+04  4.71e+04  1.31e-01  1.09e-01  5.47e-01
1.31e+05  1.01e+05  3.73e+04  1.20e-01  1.01e-01  5.47e-01
1.26e+05  1.00e+05  3.18e+04  9.77e-02  9.09e-02  5.47e-01
1.15e+05  9.11e+04  2.93e+04  7.10e-02  6.81e-02  5.47e-01
1.11e+05  8.86e+04  2.77e+04  5.80e-02  5.52e-02  5.47e-01
1.13e+05  9.20e+04  2.67e+04  5.86e-02  4.44e-02  5.47e-01
1.11e+05  9.01e+04  2.57e+04  5.25e-02  3.74e-02  6.28e-01
1.05e+05  8.59e+04  2.43e+04  3.99e-02  4.06e-02  7.44e-01
1.02e+05  8.39e+04  2.30e+04  3.24e-02  3.73e-02  7.44e-01
1.04e+05  8.54e+04  2.28e+04  3.34e-02  2.62e-02  7.44e-01
1.04e+05  8.56e+04  2.26e+04  3.23e-02  2.12e-02  8.40e-01
1.02e+05  8.42e+04  2.21e+04  2.71e-02  2.73e-02  1.04e+00
1.00e+05  8.29e+04  2.18e+04  2.31e-02  2.71e-02  1.04e+00
1.01e+05  8.32e+04  2.21e+04  2.27e-02  2.04e-02  1.04e+00
1.01e+05  8.35e+04  2.23e+04  2.26e-02  1.57e-02  1.04e+00
1.01e+05  8.29e+04  2.21e+04  2.04e-02  1.83e-02  1.24e+00
9.98e+04  8.24e+04  2.18e+04  1.83e-02  1.91e-02  1.24e+00
9.98e+04  8.24e+04  2.18e+04  1.75e-02  1.60e-02  1.24e+00
9.99e+04  8.25e+04  2.18e+04  1.72e-02  1.29e-02  1.24e+00
9.97e+04  8.24e+04  2.17e+04  1.62e-02  1.36e-02  1.43e+00
9.95e+04  8.21e+04  2.17e+04  1.51e-02  1.40e-02  1.43e+00
9.94e+04  8.21e+04  2.17e+04  1.44e-02  1.26e-02  1.43e+00
9.95e+04  8.21e+04  2.17e+04  1.40e-02  1.08e-02  1.43e+00
9.94e+04  8.20e+04  2.17e+04  1.34e-02  1.09e-02  1.63e+00
9.92e+04  8.19e+04  2.17e+04  1.27e-02  1.17e-02  1.80e+00
9.92e+04  8.18e+04  2.17e+04  1.21e-02  1.13e-02  1.80e+00
9.91e+04  8.18e+04  2.17e+04  1.18e-02  1.02e-02  1.80e+00
9.91e+04  8.18e+04  2.17e+04  1.14e-02  9.40e-03  1.80e+00
9.91e+04  8.17e+04  2.17e+04  1.09e-02  9.62e-03  1.98e+00
9.90e+04  8.17e+04  2.17e+04  1.05e-02  9.41e-03  1.98e+00
9.90e+04  8.17e+04  2.17e+04  1.01e-02  8.80e-03  1.98e+00
9.90e+04  8.16e+04  2.17e+04  9.79e-03  8.26e-03  1.98e+00
9.89e+04  8.16e+04  2.17e+04  9.46e-03  7.88e-03  1.98e+00
9.89e+04  8.16e+04  2.17e+04  9.09e-03  7.99e-03  2.17e+00
9.89e+04  8.16e+04  2.17e+04  8.78e-03  7.79e-03  2.17e+00
9.89e+04  8.15e+04  2.17e+04  8.49e-03  7.38e-03  2.17e+00
9.88e+04  8.15e+04  2.17e+04  8.20e-03  7.02e-03  2.17e+00
9.88e+04  8.15e+04  2.17e+04  7.93e-03  6.72e-03  2.17e+00
9.88e+04  8.15e+04  2.17e+04  7.67e-03  6.44e-03  2.17e+00
9.88e+04  8.15e+04  2.16e+04  7.41e-03  6.16e-03  2.17e+00
9.88e+04  8.14e+04  2.16e+04  7.14e-03  6.24e-03  2.39e+00
9.87e+04  8.14e+04  2.16e+04  6.89e-03  6.12e-03  2.39e+00
9.87e+04  8.14e+04  2.16e+04  6.66e-03  5.90e-03  2.39e+00
9.87e+04  8.14e+04  2.16e+04  6.44e-03  5.66e-03  2.39e+00
9.87e+04  8.14e+04  2.16e+04  6.22e-03  5.42e-03  2.39e+00
9.87e+04  8.14e+04  2.16e+04  6.01e-03  5.19e-03  2.39e+00
9.87e+04  8.14e+04  2.16e+04  5.81e-03  4.97e-03  2.39e+00
9.87e+04  8.14e+04  2.16e+04  5.62e-03  4.77e-03  2.39e+00
9.86e+04  8.13e+04  2.16e+04  5.43e-03  4.60e-03  2.39e+00
9.86e+04  8.13e+04  2.16e+04  5.26e-03  4.42e-03  2.39e+00
9.86e+04  8.13e+04  2.16e+04  5.09e-03  4.26e-03  2.39e+00
9.86e+04  8.13e+04  2.16e+04  4.93e-03  4.10e-03  2.39e+00
----------------------------------------------------------------

Display solve time and denoising performance.

print("TVL1Denoise solve time: %5.2f s" % b.timer.elapsed('solve'))
print("Noisy image PSNR:    %5.2f dB" % metric.psnr(img, imgn))
print("Denoised image PSNR: %5.2f dB" % metric.psnr(img, imgr))

TVL1Denoise solve time:  1.24 s
Noisy image PSNR:    12.02 dB
Denoised image PSNR: 29.29 dB

Display reference, corrupted, and denoised images.

fig = plot.figure(figsize=(20, 5))
plot.subplot(1, 3, 1)
plot.imview(img, title='Reference', fig=fig)
plot.subplot(1, 3, 2)
plot.imview(imgn, title='Corrupted', fig=fig)
plot.subplot(1, 3, 3)
plot.imview(imgr, title=r'Restored ($\ell_1$-TV)', 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()
ObjFun = [float(x) for x in its.ObjFun]
PrimalRsdl = [float(x) for x in its.PrimalRsdl]
DualRsdl = [float(x) for x in its.DualRsdl]
fig = plot.figure(figsize=(20, 5))
plot.subplot(1, 3, 1)
plot.plot(ObjFun, xlbl='Iterations', ylbl='Functional', fig=fig)
plot.subplot(1, 3, 2)
plot.plot(np.vstack((PrimalRsdl, 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()