Numerical Tours
← All Python tours
Download notebook Open in Colab

Volumetric Wavelet Processing

Extend multiscale image processing to a three-dimensional volume. Examine slices and isosurfaces, then study how coefficient thresholding changes both the interior signal and the geometry of its level sets.

Run this tour

Run the cells in order with a Python 3 kernel. The first cell locates the companion data and toolbox and installs missing dependencies when needed. All worked examples include their implementation directly in this notebook. Random seeds make comparisons reproducible; you can change them to explore other samples.

# Locate the companion toolbox locally, or fetch it for a standalone/Colab copy.
from pathlib import Path
import importlib.util
import os
import subprocess
import sys

working = Path.cwd()
candidates = [working, working / "python", working.parent / "python"]
python_dir = next((p for p in candidates if (p / "nt_toolbox").is_dir()), None)
if python_dir is None:
    checkout = working / "numerical-tours-support"
    if not checkout.exists():
        subprocess.run(
            [
                "git",
                "clone",
                "--depth",
                "1",
                "--branch",
                "master",
                "https://github.com/gpeyre/numerical-tours.git",
                str(checkout),
            ],
            check=True,
        )
    python_dir = checkout / "python"
os.chdir(python_dir)
if str(python_dir) not in sys.path:
    sys.path.insert(0, str(python_dir))
requirements = python_dir / "requirements.txt"
if any(
    importlib.util.find_spec(name) is None
    for name in [
        "numpy",
        "scipy",
        "matplotlib",
        "skimage",
        "sklearn",
        "pywt",
        "ipywidgets",
        "cvxpy",
        "skfmm",
        "autograd",
        "progressbar",
        "celer",
    ]
):
    subprocess.run(
        [sys.executable, "-m", "pip", "install", "-r", str(requirements)], check=True
    )

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(0)
plt.rcParams.update(
    {
        "figure.figsize": (8, 4),
        "figure.dpi": 100,
        "axes.spines.top": False,
        "axes.spines.right": False,
        "font.size": 11,
        "image.cmap": "gray",
    }
)
%matplotlib inline

$\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\enscond}[2]{\lbrace #1, #2 \rbrace}$ $\newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\umax}[1]{\underset{#1}{\max}\;}$ $\newcommand{\uargmin}[1]{\underset{#1}{argmin}\;}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\diag}[1]{{diag}\left( #1 \right)}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\qarrq}{ \quad \Longrightarrow \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\CC}{\mathbb{C}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\Zz}{\mathcal{Z}}$ $\newcommand{\Ww}{\mathcal{W}}$ $\newcommand{\Vv}{\mathcal{V}}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\NN}{\mathcal{N}}$ $\newcommand{\Hh}{\mathcal{H}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\Ee}{\mathcal{E}}$ $\newcommand{\Cc}{\mathcal{C}}$ $\newcommand{\Gg}{\mathcal{G}}$ $\newcommand{\Ss}{\mathcal{S}}$ $\newcommand{\Pp}{\mathcal{P}}$ $\newcommand{\Ff}{\mathcal{F}}$ $\newcommand{\Xx}{\mathcal{X}}$ $\newcommand{\Mm}{\mathcal{M}}$ $\newcommand{\Ii}{\mathcal{I}}$ $\newcommand{\Dd}{\mathcal{D}}$ $\newcommand{\Ll}{\mathcal{L}}$ $\newcommand{\Tt}{\mathcal{T}}$ $\newcommand{\si}{\sigma}$ $\newcommand{\al}{\alpha}$ $\newcommand{\la}{\lambda}$ $\newcommand{\ga}{\gamma}$ $\newcommand{\Ga}{\Gamma}$ $\newcommand{\La}{\Lambda}$ $\newcommand{\Si}{\Sigma}$ $\newcommand{\be}{\beta}$ $\newcommand{\de}{\delta}$ $\newcommand{\De}{\Delta}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$

This numerical tour explores volumetric (3D) data processing.

import numpy as np
import scipy as scp
import pylab as pyl
import matplotlib.pyplot as plt

from nt_toolbox.general import circshift, clamp, np, plt, pylab, rescale
from nt_toolbox.signal import imageplot, np, plt, pylab, snr

import warnings

%matplotlib inline

3D Volumetric Datasets

We load a volumetric data.

from nt_toolbox.read_bin import np, read_bin

M = read_bin("nt_toolbox/data/vessels.bin", ndims=3)

M = rescale(M)

Size of the image (here it is a cube).

n = np.shape(M)[1]

We can display some horizontal slices.

slices = np.round(np.linspace(10, n - 10, 4))
plt.figure(figsize=(10, 10))

for i in range(len(slices)):
    s = int(slices[i])
    imageplot(M[:, :, s], "Z = %i" % s, [2, 2, i + 1])
No description has been provided for this image

We can display an isosurface of the dataset (here we sub-sample to speed up the computation).

from nt_toolbox.isosurface import isosurface, np, plt

isosurface(M, 0.5, 3)
No description has been provided for this image

3D Haar Transform

An isotropic 3D Haar transform recursively extracts details wavelet coefficients by performing local averages/differences along the X/Y/Z axis.

We apply a step of Haar transform in the X/Y/Z direction

Initialize the transform

MW = np.copy(M)

Average/difference along X

MW = np.concatenate(
    (
        (MW[0:n:2, :, :] + MW[1:n:2, :, :]) / np.sqrt(2),
        (MW[0:n:2, :, :] - MW[1:n:2, :, :]) / np.sqrt(2),
    ),
    0,
)

Average/difference along Y

MW = np.concatenate(
    (
        (MW[:, 0:n:2, :] + MW[:, 1:n:2, :]) / np.sqrt(2),
        (MW[:, 0:n:2, :] - MW[:, 1:n:2, :]) / np.sqrt(2),
    ),
    1,
)

Average/difference along Z

MW = np.concatenate(
    (
        (MW[:, :, 0:n:2] + MW[:, :, 1:n:2]) / np.sqrt(2),
        (MW[:, :, 0:n:2] - MW[:, :, 1:n:2]) / np.sqrt(2),
    ),
    2,
)

Display a horizontal and vertical slice to see the structure of the coefficients.

plt.figure(figsize=(10, 5))
imageplot(MW[:, :, 30], "Horizontal slice", [1, 2, 1])
imageplot((MW[:, 30, :]), "Vertical slice", [1, 2, 2])
No description has been provided for this image

Worked example 1

Implement the forward wavelet transform by iteratively applying these transform steps to the low pass residual.

MW = np.copy(M)
for j in range(1, int(np.log2(n)) + 1):
    p = int(n / 2 ** (j - 1))
    sel = np.arange(0, p)
    even = np.arange(0, p, 2)
    odd = np.arange(1, p, 2)
    # average/ difference along X
    MW[np.ix_(sel, sel, sel)] = np.concatenate(
        (
            (MW[np.ix_(even, sel, sel)] + MW[np.ix_(odd, sel, sel)]) / np.sqrt(2),
            (MW[np.ix_(even, sel, sel)] - MW[np.ix_(odd, sel, sel)]) / np.sqrt(2),
        ),
        0,
    )
    # average/ difference along Y
    MW[np.ix_(sel, sel, sel)] = np.concatenate(
        (
            (MW[np.ix_(sel, even, sel)] + MW[np.ix_(sel, odd, sel)]) / np.sqrt(2),
            (MW[np.ix_(sel, even, sel)] - MW[np.ix_(sel, odd, sel)]) / np.sqrt(2),
        ),
        1,
    )
    # average/ difference along Z
    MW[np.ix_(sel, sel, sel)] = np.concatenate(
        (
            (MW[np.ix_(sel, sel, even)] + MW[np.ix_(sel, sel, odd)]) / np.sqrt(2),
            (MW[np.ix_(sel, sel, even)] - MW[np.ix_(sel, sel, odd)]) / np.sqrt(2),
        ),
        2,
    )

Volumetric Data Haar Approximation

An approximation is obtained by keeping only the largest coefficients.

We threshold the coefficients to perform $m$-term approximation.

number of kept coefficients

from nt_toolbox.perform_thresholding import np, perform_thresholding

m = round(0.01 * n**3)
MWT = perform_thresholding(MW, m, type="largest")

Worked example 2

Implement the backward transform to compute an approximation $M_1$ from the coefficients MWT.

M1 = np.copy(MWT)

for j in range(int(np.log2(n)), 0, -1):
    p = int(n / 2**j)
    sel = np.arange(0, p)
    sel1 = np.arange(0, 2 * p)
    selw = np.arange(p, 2 * p)
    even = np.arange(0, 2 * p, 2)
    odd = np.arange(1, 2 * p, 2)
    # average/ difference along X
    A = M1[np.ix_(sel, sel1, sel1)]
    D = M1[np.ix_(selw, sel1, sel1)]
    M1[np.ix_(even, sel1, sel1)] = (A + D) / np.sqrt(2)
    M1[np.ix_(odd, sel1, sel1)] = (A - D) / np.sqrt(2)
    # average/ difference along Y
    A = M1[np.ix_(sel1, sel, sel1)]
    D = M1[np.ix_(sel1, selw, sel1)]
    M1[np.ix_(sel1, even, sel1)] = (A + D) / np.sqrt(2)
    M1[np.ix_(sel1, odd, sel1)] = (A - D) / np.sqrt(2)
    # average/ difference along Z
    A = M1[np.ix_(sel1, sel1, sel)]
    D = M1[np.ix_(sel1, sel1, selw)]
    M1[np.ix_(sel1, sel1, even)] = (A + D) / np.sqrt(2)
    M1[np.ix_(sel1, sel1, odd)] = (A - D) / np.sqrt(2)

Display the approximation as slices.

s = 30

plt.figure(figsize=(10, 5))
imageplot(M[:, :, s], "Original", [1, 2, 1])
imageplot(clamp(M1[:, :, s]), "Approximation", [1, 2, 2])
No description has been provided for this image

Display the approximated isosurface.

isosurface(M1, 0.5, 2)
No description has been provided for this image

Linear Volumetric Denoising

Linear denoising is obtained by low pass filtering.

We add a Gaussian noise to the image.

from numpy import random

sigma = 0.06
Mnoisy = M + sigma * random.randn(n, n, n)

Display slices of the noisy data.

plt.figure(figsize=(10, 5))
imageplot(Mnoisy[:, :, n // 2], "X slice", [1, 2, 1])
imageplot(Mnoisy[:, n // 2, :], "Y slice", [1, 2, 2])
No description has been provided for this image

A simple denoising method performs a linear filtering of the data.

We build a Gaussian filter of width $\sigma$.

Construct a 3D grid

x = np.arange(-n // 2, n // 2)
[X, Y, Z] = np.meshgrid(x, x, x)

Gaussian filter

s = 2  # width
h = np.exp(-(X**2 + Y**2 + Z**2) / (2 * s**2))
h = h / np.sum(h)

The filtering is computed over the Fourier domain.

Mh = np.real(
    pyl.ifft2(
        pyl.fft2(Mnoisy, axes=(0, 1, 2))
        * pyl.fft2(pyl.fftshift(h, axes=(0, 1, 2)), axes=(0, 1, 2)),
        axes=(0, 1, 2),
    )
)

Display denoised slices.

i = 40
plt.figure(figsize=(10, 5))
imageplot(Mnoisy[:, :, i], "Noisy", [1, 2, 1])
imageplot(Mh[:, :, i], "Denoised", [1, 2, 2])
No description has been provided for this image

Display denoised iso-surface.

isosurface(M, 0.5, 3)
No description has been provided for this image

Worked example 3

Select the optimal blurring width $s$ to reach the smallest possible SNR. Keep the optimal denoising Mblur.

ntests = 20
slist = np.linspace(0.01, 1.5, ntests)
err = []
for i in range(ntests):
    h = np.exp(-(X**2 + Y**2 + Z**2) / (2 * slist[i] ** 2))
    h = h / np.sum(h)
    Mh = np.real(
        pyl.ifft2(
            pyl.fft2(Mnoisy, axes=(0, 1, 2))
            * pyl.fft2(pyl.fftshift(h, axes=(0, 1, 2)), axes=(0, 1, 2)),
            axes=(0, 1, 2),
        )
    )
    err = err + [snr(M, Mh)]
    if i > 1 and err[i] > np.max(err[:i]):
        Mblur = Mh

plt.figure(figsize=(7, 5))
plt.plot(slist, err, ".-")
plt.xlabel("s")
plt.ylabel("SNR")
plt.show()
No description has been provided for this image

Display optimally denoised iso-surface.

isosurface(Mblur, 0.5, 2, "Filtering, SNR = %.1f dB" % snr(M, Mblur))
No description has been provided for this image

Non-Linear Wavelet Volumetric Denoising

Denoising is obtained by removing small amplitude coefficients that corresponds to noise.

Worked example 4

Perforn Wavelet denoising by thresholding the wavelet coefficients of Mnoisy. Test both hard thresholding and soft thresholding to determine the optimal threshold and the corresponding SNR. Record the optimal result Mwav.

from nt_toolbox.perform_haar_transf import np, perform_haar_transf
from nt_toolbox.perform_thresholding import np, perform_thresholding

MW = perform_haar_transf(Mnoisy, 1, +1)
Tlist = np.linspace(1, 4, 20) * sigma
err_hard = []
err_soft = []

for i in range(len(Tlist)):
    MWT = perform_thresholding(MW, Tlist[i], "hard")
    M1 = perform_haar_transf(MWT, 1, -1)
    err_hard = err_hard + [snr(M, M1)]
    MWT = perform_thresholding(MW, Tlist[i], "soft")
    M1 = perform_haar_transf(MWT, 1, -1)
    err_soft = err_soft + [snr(M, M1)]
    if i > 1 and err_soft[i] > np.max(err_soft[:i]):
        Mwav = M1

plt.figure(figsize=(7, 5))

plt.plot(Tlist / sigma, err_hard, ".-", label="hard", color="b")
plt.plot(Tlist / sigma, err_soft, ".-", label="soft", color="r")
plt.xlabel(r"$T/\sigma$")
plt.ylabel("SNR")
plt.ylim(np.min(err_hard), np.max(err_soft))
plt.legend(loc="upper right")

plt.show()
No description has been provided for this image

Display denoised iso-surface with optimal soft thresholding.

isosurface(Mblur, 0.5, 2, "Soft thresholding, SNR = %.1f dB" % snr(M, Mwav))
No description has been provided for this image

Orthogonal wavelet thresholdings suffers from blocking artifacts. This can be aleviated by performing a cycle spinning denoising, which averages the denosing result of translated version of the signal.

A typical cycle spinning process is like this.

Maximum translation.

w = 4

List of translations.

[dZ, dX, dY] = np.meshgrid(np.arange(0, w), np.arange(0, w), np.arange(0, w))
dX = np.ravel(dX)
dY = np.ravel(dY)
dZ = np.ravel(dZ)

Initialize spinning process.

Mspin = np.zeros([n, n, n])

Spin.

def circshift(x, v):
    x = np.roll(x, v[0], axis=0)
    x = np.roll(x, v[1], axis=1)
    x = np.roll(x, v[2], axis=2)
    return x


for i in range(w**3):
    # shift the image
    MnoisyC = circshift(Mnoisy, [dX[i], dY[i], dZ[i]])
    # denoise the image to get a result M1
    M1 = MnoisyC  # replace this line by some denoising
    # shift inverse
    M1 = circshift(M1, [-dX[i], -dY[i], -dZ[i]])
    # average the result
    Mspin = Mspin * (i) / (i + 1) + M1 / (i + 1)

Worked example 5

Implement cycle spinning hard thresholding with $T=3\sigma$.

T = 3 * sigma
w = 4
Mspin = np.zeros([n, n, n])

for i in range(w**3):
    # shift the image
    MnoisyC = circshift(Mnoisy, [dX[i], dY[i], dZ[i]])
    # denoise
    MW = perform_haar_transf(MnoisyC, 1, +1)
    MWT = perform_thresholding(MW, T, "hard")
    M1 = perform_haar_transf(MWT, 1, -1)
    # back
    M1 = circshift(M1, [-dX[i], -dY[i], -dZ[i]])
    # average the result
    Mspin = Mspin * (i) / (i + 1) + M1 / (i + 1)

Display denoised iso-surface.

isosurface(Mspin, 0.5, 2, "Cycle spinning, SNR = %.1f dB" % snr(M, Mspin))
No description has been provided for this image

References and further reading