Show code cell content Hide code cell content

import mmf_setup;mmf_setup.nbinit()
import logging;logging.getLogger('matplotlib').setLevel(logging.CRITICAL)
%matplotlib inline
import numpy as np, matplotlib.pyplot as plt

$\newcommand{\vect}[1]{\mathbf{#1}} \newcommand{\uvect}[1]{\hat{#1}} \newcommand{\abs}[1]{\lvert#1\rvert} \newcommand{\norm}[1]{\lVert#1\rVert} \newcommand{\I}{\mathrm{i}} \newcommand{\ket}[1]{\left|#1\right\rangle} \newcommand{\bra}[1]{\left\langle#1\right|} \newcommand{\braket}[1]{\langle#1\rangle} \newcommand{\Braket}[1]{\left\langle#1\right\rangle} \newcommand{\op}[1]{\mathbf{#1}} \newcommand{\mat}[1]{\mathbf{#1}} \newcommand{\d}{\mathrm{d}} \newcommand{\D}[1]{\mathcal{D}[#1]\;} \newcommand{\pdiff}[3][]{\frac{\partial^{#1} #2}{\partial {#3}^{#1}}} \newcommand{\diff}[3][]{\frac{\d^{#1} #2}{\d {#3}^{#1}}} \newcommand{\ddiff}[3][]{\frac{\delta^{#1} #2}{\delta {#3}^{#1}}} \newcommand{\floor}[1]{\left\lfloor#1\right\rfloor} \newcommand{\ceil}[1]{\left\lceil#1\right\rceil} % The following allows KaTeX to render these for significantly improved % performance on CoCalc. \newcommand{\DeclareMathOperator}[2]{\newcommand{#1}{#2}} \DeclareMathOperator{\Tr}{Tr} \DeclareMathOperator{\erf}{erf} \DeclareMathOperator{\erfi}{erfi} \DeclareMathOperator{\sech}{sech} \DeclareMathOperator{\sn}{sn} \DeclareMathOperator{\cn}{cn} \DeclareMathOperator{\dn}{dn} \DeclareMathOperator{\sgn}{sgn} \DeclareMathOperator{\order}{O} \DeclareMathOperator{\diag}{diag} % \newcommand{\mylabel}[1]{\label{#1}\tag{#1}} \newcommand{\degree}{\circ}$

This cell adds /home/docs/checkouts/readthedocs.org/user_builds/iscimath-583-learning-from-signals/checkouts/latest/src to your path, and contains some definitions for equations and some CSS for styling the notebook. If things look a bit strange, please try the following:

• Choose "Trust Notebook" from the "File" menu.
• Re-execute this cell.
• Reload the notebook.

Code

Here we describe the various tools implemented in the math_583 package.

math_583

The module provides the following classes:

Math 583 - Learning from Signals

This package contains source code for the Spring 2023 offering of Math 583: iSciMath - Learning from Images and Signals: Theory and Computation at Washington State University (WSU).

math_583.denoise

This model is a proxy for

\[\begin{multline*} E[u] = \int\frac{1}{p}\abs{\vect{\nabla}u}^p \\ + \frac{λ}{q}\int\abs{u - d}^q,\\ E'[u] = -\vect{\nabla}\cdot\Bigl(\vect{\nabla}u \abs{\vect{\nabla}u}^{p-2}\Bigr) \\ + λ(u-d)\abs{u-d}^{q-2}. \end{multline*}\]

These are formally equivalent in the continuum limit with vanishing regulators \(ϵ_p, ϵ_q \rightarrow 0\) with an appropriate derivative operator, but numerically distinct as discussed in Product rule.

This module provides tools to denoise images based on the following model:

\[\begin{gather*} E[u] = \int\frac{1}{p}\abs{\vect{\nabla}u}^p + \frac{λ}{q}\int\abs{u - d}^q. \end{gather*}\]

Minimizing is equivalent to solving:

\[\begin{gather*} E'[u] = -\nabla^2u + λ(u - d) = 0. \end{gather*}\]

A direct approach to the minimization procedure is given by Denoise.minimize() which uses scipy.optimize.minimize() and the L-BFGS-B method). Minimization is not the fastest, but is quite general, and be directly applied to more complicated methods.

\[\begin{gather*} E[u] = \int\frac{1}{p}\bigl(\abs{\vect{\nabla}u}^2 + ϵ_p\bigr)^{p/2} + \frac{λ}{q}\int\bigl((u - d)^2 + ϵ_q\bigr)^{q/2},\\ E'[u] = -\vect{\nabla}\cdot\Bigl( (\vect{\nabla} u) \bigl( \abs{\vect{\nabla}u}^2 + ϵ_p \bigr)^{(p-2)/2} \Bigr) + λ(u - d)\bigl((u - d)^2 + ϵ_q\bigr)^{(q-2)/2}\\ \end{gather*}\]

Details about anti-symmetric derivatives \(\mat{D}^T = - \mat{D}\).

In order to obtain the expression for \(E'[u]\), we need to compute

\[\begin{gather*} \pdiff{\bigl[\abs{\vect{\nabla}u}^2\bigr]_{a}}{u_b} \end{gather*}\]

which appears in the first integral after the chain rule is applied. Noting that \(\vect{\nabla}\) is a linear operator, we will have a matrix representation \(\vect{\mat{D}}_{ac}\) such that

\[\begin{gather*} \bigl[\abs{\vect{\nabla} u}^2\bigr]_{a} = \sum_{i, c, d}D^{i}_{ac}u_c D^{i}_{ad}u_d\\ \pdiff{\bigl[\abs{\vect{\nabla}u}^2\bigr]_{a}}{u_b} = \sum_{i, d}D^{i}_{ab} D^{i}_{ad}u_{d} + \sum_{i, c}D^{i}_{ac}u_{c} D^{i}_{ab} = 2\sum_{i, c}D^{i}_{ab} D^{i}_{ac}u_{c}. \end{gather*}\]

When included in the integral with an additional factor \(f_{a} = \bigl(\abs{\vect{\nabla}u}^2 + ϵ_p\bigr)^{(p-2)/2}_{a}\) we end up with

\[\begin{gather*} 2\sum_{i, a, c}D^{i}_{ab} D^{i}_{ac}u_{c} f_{a} = 2 \sum_{a} \underbrace{(-D^{i}_{ba})}_{D^{i}_{ab}}\Bigl(\sum_{i, c} D^{i}_{ac}u_{c} f_{a}\Bigr) \equiv - 2 \vect{\nabla}\cdot\Bigl((\vect{\nabla}u) f\Bigr), \end{gather*}\]

iff the matrix representations of the derivatives are anti-symmetric

\[\begin{gather*} \mat{D}^{\dagger} = - \mat{D}, \qquad D^{i}_{ab} = - D^{i}_{ba}. \end{gather*}\]

Alternatively, the outer derivative must be computed using the transpose of the original derivative matrix as in the explicit formula. Note: In the above, we have not consider complex values. Similar formulae apply with appropriate complex conjugations applied.

A faster approach is provided by Denoise.solve(), which directly solves the minimization condition using Fourier Techniques.

Tools for exploring image denoising.

subplots(cols=1, rows=1, height=3, aspect=1, **kw)

More convenient subplots that automatically sets the figsize.

Parameters

• cols (int) – Number of columns and rows in the figure.

• rows (int) – Number of columns and rows in the figure.

• height (float) – Height of an individual element. Unless specified, the figure size will be figsize=(cols * height * aspect, rows * height).

• aspect (float) – Aspect ratio of each figure (width/height).

• **kw (dict) – Other arguments are passed to plt.subplot()

class Image(data=None, **kw)

Class to load and process images.

data

If provided, then this is used as the image data. If the first dimension has length 3 or 4, then it is assumed to be RGB or RGBA data respectively. If the dtype is np.uint8, then it is unconverted, otherwise, it is normalized using the default colormap.

Type

array_like | None

dir

Directory with images.

Type

Path

filename

Location of file to be loaded if data is None.

Type

bytes

seed

Seed for random number generator.

Type

int

property rgb

Return the RGB form of the image.

get_data(normalize=True, sigma=0, rng=None)

Return greyscale image.

Parameters

• normalize (bool) – If True, then normalize the data so it is between 0 and 1.

• sigma (float) – Standard deviation (as a fraction) of gaussian noise to add to the image. The result will be clipped so it does not exceed (0, 255) or (0, 1) if normalize==True.

class Denoise(image=None, **kw)

Class for denoising images.

image

Instance of Image with the image data.

Type

Image

lam

Parameter λ controlling the regularization. Larger values will produce an image closer to the target. Smaller values will produce smoother images.

Type

float

p, q

Powers in the regularization term and the data fidelity terms respectively.

Type

float

use_shortcuts

If True, then use specialized shortcuts where applicable, otherwise use the general code.

Type

bool

eps_p, eps_q

Regularization constants for the regularization term and the data fidelity terms respectively.

Type

float

real

If True, then some operations that might be complex will return real values.

Type

bool

laplacian(u)

Return the laplacian of u.

derivative1d(input, axis=- 1, output=None, mode=None, cval=0.0, kind='forward', compute=False, transpose=False)

Return the difference of input along the specified axis.

This is intended to be used as the derivative function in various scipy.ndimage routines.

Parameters

• kind ('forward', 'backward', 'centered') –

• transpose (Bool) – If True, then apply the negative transpose of the derivative operator.

gradient(u, real=True)

Return the gradient of u.

Parameters

real (bool) – If True, then return only the real part. It is important to set this to False if computing the Laplacian as divergence(gradient(u)) for example. Only applies for mode == “periodic”.

divergence(v, transpose=True)

Return the divergence of v.

gradient_magnitude(u)

Return the absolute magnitude of the gradient of u.

gradient_magnitude2(u)

Return the square of the magnitude of the gradient of u.

get_energy(u, parts=False, normalize=False)

Return the energy.

Parameters

• parts (bool) – If True, return (E, E_regularization, E_data_fidelity)

• normalize (bool) – If True, normalize by the starting values for u_noise.

get_denergy(u, normalize=False)

Return E’(u).

pack(u)

Return y, the 1d real representation of u for solving.

unpack(y)

Return u from the 1d real representation y.

compute_dy_dt(t, y)

Return dy_dt for the solver.

minimize(u0=None, method='L-BFGS-B', callback=True, tol=1e-08, plot=False, **kw)

Directly solve the minimization problem with the L-BFGS-B method.

solve()

Directly solve the minimization problem using Fourier techniques.

This is much faster than running minimize but only supports limited modes.