polsar-h-alpha-plane.py

"""
S. R. Cloude and E. Pottier, "An entropy based classification scheme for land applications of polarimetric SAR," in IEEE Transactions on Geoscience and Remote Sensing, vol. 35, no. 1, pp. 68-78, Jan. 1997, doi: 10.1109/36.551935.
"""
import matplotlib.pyplot as plt
import numpy as np


def h_alpha_decomposition(T3):
    assert isinstance(T3, np.ndarray)
    assert T3.ndim >= 2
    assert T3.shape[-2:] == (3, 3)

    EPS = 1e-30
    eig_values, eig_vectors = np.linalg.eigh(T3)
    eig_values[eig_values < EPS] = EPS
    probs = eig_values / np.sum(eig_values, axis=-1, keepdims=True)
    h = -np.sum(probs * (np.log(probs) / np.log(3)), axis=-1)
    alpha = np.sum(probs * np.degrees(np.arccos(np.abs(eig_vectors[..., 0, :]))), axis=-1)
    return h, alpha


def T3_curve1(m):
    assert 0.0 <= m <= 1.0
    return np.array(
        [
            [1, 0, 0],
            [0, m, 0],
            [0, 0, m],
        ]
    )


def T3_curve2(m):
    assert 0.0 <= m <= 1.0
    return (
        np.array(
            [
                [0, 0, 0],
                [0, 1, 0],
                [0, 0, 2 * m],
            ]
        )
        if m < 0.5
        else np.array(
            [
                [2 * m - 1, 0, 0],
                [0, 1, 0],
                [0, 0, 1],
            ]
        )
    )


if __name__ == "__main__":
    fig, ax = plt.subplots(figsize=(5, 5))

    curve_style = {
        "color": "black",
        "linewidth": 1,
    }
    ax.plot(*h_alpha_decomposition(np.array([T3_curve1(m) for m in np.linspace(0, 1, 100)])), **curve_style)
    ax.plot(*h_alpha_decomposition(np.array([T3_curve2(m) for m in np.linspace(0, 1, 100)])), **curve_style)

    bounds = [
        ([0.0, 0.5], [42.5, 42.5]),
        ([0.0, 0.5], [47.5, 47.5]),
        ([0.5, 0.5], [0.0, 90.0]),
        ([0.5, 0.9], [40.0, 40.0]),
        ([0.5, 0.9], [50.0, 50.0]),
        ([0.9, 0.9], [0.0, 90.0]),
        ([0.9, 1.0], [40.0, 40.0]),
        ([0.9, 1.0], [55.0, 55.0]),
    ]
    for xs, ys in bounds:
        ax.plot(xs, ys, "--", color="gray", linewidth=1)

    ax.set_xlim(0, 1)
    ax.set_ylim(0, 90)
    ax.set_xlabel("Entropy")
    ax.set_ylabel("Alpha")
    ax.tick_params(top="on", right="on", direction="in")
    plt.show()