This page was generated from docs/notebooks/others/derivative_operator.ipynb.
Download notebook.

Derivative Operation¶

Here we show you examples of derivative/laplacian matrices and how they work on the example image.

Please see the derivative theory page if you want to know the details of the theory.

[1]:

import numpy as np
import ultraplot as uplt
from matplotlib import pyplot as plt
from matplotlib.cbook import get_sample_data
from PIL import Image

from cherab.inversion.tools import derivative_matrix, laplacian_matrix

Visualize derivative matrix¶

Let us compute the simple derivative matrix applying a 3 x 3 image array.

Firstly we show derivative matrices along to vertical direction (axis=0 direction) with different schemes. The derivative matrices are shown by plotting the matrix elements in 2-D image.

[2]:

fig, axes = uplt.subplots(ncols=3)
for ax, scheme in zip(axes, ["forward", "backward", "central"], strict=True):
    # compute derivative matrix
    dmat = derivative_matrix((3, 3), axis=0, scheme=scheme).toarray()
    dmat[dmat == 0] = np.nan

    # plot derivative matrix
    ax.heatmap(
        dmat,
        cmap="ColdHot",
        vmin=-1,
        vmax=1,
        lw=0.1,
        ec="k",
        labels=True,
        precision=2,
        clip_on=False,
    )
    ax.format(title=scheme)

axes.format(
    xloc="top",
    yloc="left",
    yreverse=True,
    ticklabelweight="bold",
    linewidth=0,
    tickpad=4,
)

../../_images/notebooks_others_derivative_operator_5_0.png

Also we show the horizontal derivative matrix (axis=1 direction) case.

[3]:

fig, axes = uplt.subplots(ncols=3)
for ax, scheme in zip(axes, ["forward", "backward", "central"], strict=True):
    # compute derivative matrix
    dmat = derivative_matrix((3, 3), axis=1, scheme=scheme).toarray()
    dmat[dmat == 0] = np.nan

    # plot derivative matrix
    ax.heatmap(
        dmat,
        cmap="ColdHot",
        vmin=-1,
        vmax=1,
        lw=0.1,
        ec="k",
        labels=True,
        precision=2,
        clip_on=False,
    )
    ax.format(title=scheme)

axes.format(
    xloc="top",
    yloc="left",
    yreverse=True,
    ticklabelweight="bold",
    linewidth=0,
    tickpad=4,
)

../../_images/notebooks_others_derivative_operator_7_0.png

We can see that some elements are zero due to the boundary condition.

Visualize laplacian matrix¶

Let us compute the simple laplacian matrix applying a 3 x 3 image array and show it by printing the matrix elements and plotting the matrix elements in 2-D image. There are two types of laplacian matrix: difference between only orthogonal neighbors or all neighbors including diagonal neighbors. The step size along diagonal direction is \(\sqrt{2}h\).

[4]:

fig, axes = uplt.subplots(ncols=2)
for ax, diagonal in zip(axes, [False, True], strict=True):
    # compute laplacian matrix
    lmat = laplacian_matrix((3, 3), diagonal=diagonal).toarray()
    lmat[lmat == 0] = np.nan

    # plot laplacian matrix
    ax.heatmap(
        lmat,
        cmap="ColdHot",
        vmin=-6,
        vmax=1,
        lw=0.1,
        ec="k",
        labels=True,
        precision=2,
        clip_on=False,
    )
    title = "diagonal" if diagonal else "no diagonal"
    ax.format(title=title)

axes.format(
    xloc="top",
    yloc="left",
    yreverse=True,
    ticklabelweight="bold",
    linewidth=0,
    tickpad=4,
)

../../_images/notebooks_others_derivative_operator_10_0.png

Apply the derivative/laplacian matrix to a sample image¶

Next, let us to apply derivative/laplacian matrices to pixels of a sample image. We use sample image stored in matplotlib library.

Let us preview the sample image converted to monochrome image.

[5]:

with get_sample_data("grace_hopper.jpg") as file:
    arr_image = plt.imread(file)

# resize the image to 1/4 of the original size
with Image.fromarray(arr_image, mode="RGB") as im:
    (width, height) = (im.width // 4, im.height // 4)
    arr_image = np.array(im.resize((width, height)))

# convert RGB image to monotonic one
arr_image = arr_image.mean(axis=2)

# show image
print(f"image array shape: {arr_image.shape}")
fig, ax = uplt.subplots()
ax.imshow(arr_image, cmap="gray")
ax.format(
    xlabel="$x$",
    ylabel="$y$",
)

image array shape: (150, 128)

../../_images/notebooks_others_derivative_operator_13_1.png

Then create the derivative/laplacian matrices. Note that the \(x\), \(y\) directions correspond to the 1, 0 axis of the image array, respectively, and the positive direction of the \(y\) axis is downward.

[6]:

# compute each derivative matrix along each axis with forward scheme
dmat_x = derivative_matrix(arr_image.shape, axis=1, scheme="forward")
dmat_y = derivative_matrix(arr_image.shape, axis=0, scheme="forward")

# compute derivative matrices along y direction with different schemes
dmat_y_backward = derivative_matrix(arr_image.shape, axis=0, scheme="backward")
dmat_y_central = derivative_matrix(arr_image.shape, axis=0, scheme="central")

# compute laplacian matrices with and without diagonal
lmat_diag = laplacian_matrix(arr_image.shape, diagonal=True)
lmat_no_diag = laplacian_matrix(arr_image.shape, diagonal=False)

Check the difference effect of the derivative matrix between \(x\) and \(y\) directions¶

[7]:

arr_derivatives = []
for dmat in [dmat_x, dmat_y]:
    arr_derivative = dmat @ arr_image.flatten()
    arr_derivatives.append(arr_derivative.reshape(arr_image.shape))

fig, axes = uplt.subplots(ncols=2, spanx=False)
for ax, arr, title in zip(axes, arr_derivatives, ["$x$", "$y$"], strict=True):
    im = ax.imshow(arr, colorbar="b")
    ax.format(title=f"{title} direction")

axes.format(
    xlabel="$x$",
    ylabel="$y$",
)

../../_images/notebooks_others_derivative_operator_17_0.png

We can see that each derivative image varies depending on each image’s gradient direction.

Check the difference derivative scheme (forward, backward, central)¶

Let us check the difference between derivative schemes in the y direction.

[8]:

arr_derivatives = []
for dmat in [dmat_y, dmat_y_backward, dmat_y_central]:
    arr_derivative = dmat @ arr_image.flatten()
    arr_derivatives.append(arr_derivative.reshape(arr_image.shape))

fig, axes = uplt.subplots(ncols=3, spanx=False)
for ax, arr, title in zip(
    axes,
    arr_derivatives,
    ["forward", "backward", "central"],
    strict=True,
):
    im = ax.imshow(arr, colorbar="b")
    ax.format(title=f"{title} direction")

axes.format(
    xlabel="$x$",
    ylabel="$y$",
)

../../_images/notebooks_others_derivative_operator_20_0.png

There is no apparent difference between the derivative schemes. The difference seems be more apparent when the gradient is large or the area near the boundary because each scheme has different boundary condition.

Check the difference effect of the laplacian matrix w/ or w/o diagonal neighbors¶

[9]:

fig, axes = uplt.subplots(ncols=2, spanx=False)
for ax, lmat, title in zip(
    axes,
    [lmat_no_diag, lmat_diag],
    ["w/o diagonal", "w/ diagonal"],
    strict=True,
):
    # compute laplacian of image
    arr_laplacian = lmat @ arr_image.flatten()
    arr_laplacian = arr_laplacian.reshape(arr_image.shape)

    # plot laplacian of image
    im = ax.imshow(
        arr_laplacian,
        colorbar="b",
    )
    ax.format(title=title)

axes.format(
    xlabel="$x$",
    ylabel="$y$",
)

../../_images/notebooks_others_derivative_operator_23_0.png

Both laplacian matrices can detect the edge of the image. The laplacian with diagonal neighbors generates larger values than that without diagonal neighbors.

In image processing, the both derivative and laplacian matrices are used for edge detection. In the inversion problem, specific to the tomography, these matrices make a reconstructed image smoother by adding the derivative/laplacian regularization term to the objective function.