.ipynb

Disk Harmonic Analysis#

import numpy as np
import pandas as pd

import matplotlib.pyplot as plt
import plotly.graph_objects as go
import seaborn as sns

from sklearn.decomposition import PCA

from ktch.datasets import load_surface_leaf_bending
from ktch.harmonic import DiskHarmonicAnalysis, xy2polar

Load the surface leaf bending dataset#

This dataset contains 60 synthetic 3D leaf surfaces with bending deformation, each paired with a disk parameterization mesh (unit disk, z=0).

data = load_surface_leaf_bending()

print(f"Number of specimens: {len(data.vertices)}")
print(f"Specimen 1 — vertices: {data.vertices[0].shape}, faces: {data.faces[0].shape}")
print(f"Specimen 1 — param vertices: {data.param_vertices[0].shape}")

Downloading data from 'https://pub-c1d6dba6c94843f88f0fd096d19c0831.r2.dev/datasets/surface_leaf_bending/v1/surface_leaf_bending.zip' to file '/home/runner/.cache/ktch-data/surface_leaf_bending/v1/surface_leaf_bending.zip'.

Number of specimens: 60
Specimen 1 — vertices: (8099, 3), faces: (15615, 3)
Specimen 1 — param vertices: (8099, 3)

df_meta = pd.DataFrame(data.meta)
df_meta["bending_angle"] = (
    np.degrees(df_meta["alphaB"]).round().astype(int).astype(str) + "\u00b0"
)
df_meta["aspect_ratio"] = df_meta["alpha"]
df_meta
lmax c alpha A alphaB kB alpha_random alphaB_random bending_angle aspect_ratio
1 0.4 0.65 0.08 1.0 2.443461 0.2 0.060572 2.528513 140° 0.08
2 0.4 0.65 0.08 1.0 2.443461 0.2 0.090911 2.318855 140° 0.08
3 0.4 0.65 0.08 1.0 2.443461 0.2 0.088957 2.417189 140° 0.08
4 0.4 0.65 0.08 1.0 2.443461 0.2 0.074596 2.506030 140° 0.08
5 0.4 0.65 0.08 1.0 2.443461 0.2 0.098990 2.294188 140° 0.08
6 0.4 0.65 0.08 1.0 2.443461 0.2 0.073196 2.492090 140° 0.08
7 0.4 0.65 0.08 1.0 2.443461 0.2 0.077120 2.287193 140° 0.08
8 0.4 0.65 0.08 1.0 2.443461 0.2 0.083971 2.280600 140° 0.08
9 0.4 0.65 0.08 1.0 2.443461 0.2 0.073029 2.431389 140° 0.08
10 0.4 0.65 0.08 1.0 2.443461 0.2 0.099243 2.553002 140° 0.08
11 0.4 0.65 0.08 1.0 1.396263 0.2 0.089302 1.409695 80° 0.08
12 0.4 0.65 0.08 1.0 1.396263 0.2 0.076330 1.242604 80° 0.08
13 0.4 0.65 0.08 1.0 1.396263 0.2 0.060501 1.299788 80° 0.08
14 0.4 0.65 0.08 1.0 1.396263 0.2 0.087227 1.536335 80° 0.08
15 0.4 0.65 0.08 1.0 1.396263 0.2 0.098797 1.548685 80° 0.08
16 0.4 0.65 0.08 1.0 1.396263 0.2 0.076723 1.280753 80° 0.08
17 0.4 0.65 0.08 1.0 1.396263 0.2 0.093675 1.558039 80° 0.08
18 0.4 0.65 0.08 1.0 1.396263 0.2 0.095911 1.566630 80° 0.08
19 0.4 0.65 0.08 1.0 1.396263 0.2 0.078236 1.387295 80° 0.08
20 0.4 0.65 0.08 1.0 1.396263 0.2 0.095689 1.447270 80° 0.08
21 0.4 0.65 0.08 1.0 0.349066 0.2 0.082873 0.390205 20° 0.08
22 0.4 0.65 0.08 1.0 0.349066 0.2 0.064031 0.266366 20° 0.08
23 0.4 0.65 0.08 1.0 0.349066 0.2 0.089630 0.373168 20° 0.08
24 0.4 0.65 0.08 1.0 0.349066 0.2 0.078860 0.330521 20° 0.08
25 0.4 0.65 0.08 1.0 0.349066 0.2 0.087673 0.318159 20° 0.08
26 0.4 0.65 0.08 1.0 0.349066 0.2 0.088925 0.305914 20° 0.08
27 0.4 0.65 0.08 1.0 0.349066 0.2 0.077160 0.372707 20° 0.08
28 0.4 0.65 0.08 1.0 0.349066 0.2 0.083355 0.268017 20° 0.08
29 0.4 0.65 0.08 1.0 0.349066 0.2 0.083693 0.406805 20° 0.08
30 0.4 0.65 0.08 1.0 0.349066 0.2 0.081946 0.311726 20° 0.08
31 0.4 0.65 0.16 1.0 2.443461 0.2 0.173997 2.326146 140° 0.16
32 0.4 0.65 0.16 1.0 2.443461 0.2 0.178168 2.367912 140° 0.16
33 0.4 0.65 0.16 1.0 2.443461 0.2 0.166143 2.497124 140° 0.16
34 0.4 0.65 0.16 1.0 2.443461 0.2 0.142324 2.521001 140° 0.16
35 0.4 0.65 0.16 1.0 2.443461 0.2 0.159544 2.492420 140° 0.16
36 0.4 0.65 0.16 1.0 2.443461 0.2 0.176739 2.296786 140° 0.16
37 0.4 0.65 0.16 1.0 2.443461 0.2 0.170455 2.304695 140° 0.16
38 0.4 0.65 0.16 1.0 2.443461 0.2 0.148710 2.478740 140° 0.16
39 0.4 0.65 0.16 1.0 2.443461 0.2 0.141903 2.555287 140° 0.16
40 0.4 0.65 0.16 1.0 2.443461 0.2 0.171313 2.610918 140° 0.16
41 0.4 0.65 0.16 1.0 1.396263 0.2 0.176521 1.244129 80° 0.16
42 0.4 0.65 0.16 1.0 1.396263 0.2 0.167533 1.252522 80° 0.16
43 0.4 0.65 0.16 1.0 1.396263 0.2 0.175097 1.438110 80° 0.16
44 0.4 0.65 0.16 1.0 1.396263 0.2 0.146855 1.570428 80° 0.16
45 0.4 0.65 0.16 1.0 1.396263 0.2 0.143005 1.422087 80° 0.16
46 0.4 0.65 0.16 1.0 1.396263 0.2 0.174511 1.497144 80° 0.16
47 0.4 0.65 0.16 1.0 1.396263 0.2 0.170496 1.470243 80° 0.16
48 0.4 0.65 0.16 1.0 1.396263 0.2 0.151212 1.282264 80° 0.16
49 0.4 0.65 0.16 1.0 1.396263 0.2 0.176779 1.257857 80° 0.16
50 0.4 0.65 0.16 1.0 1.396263 0.2 0.157964 1.505475 80° 0.16
51 0.4 0.65 0.16 1.0 0.349066 0.2 0.166466 0.346456 20° 0.16
52 0.4 0.65 0.16 1.0 0.349066 0.2 0.177635 0.386719 20° 0.16
53 0.4 0.65 0.16 1.0 0.349066 0.2 0.155776 0.441191 20° 0.16
54 0.4 0.65 0.16 1.0 0.349066 0.2 0.175787 0.407493 20° 0.16
55 0.4 0.65 0.16 1.0 0.349066 0.2 0.146865 0.504861 20° 0.16
56 0.4 0.65 0.16 1.0 0.349066 0.2 0.163440 0.445301 20° 0.16
57 0.4 0.65 0.16 1.0 0.349066 0.2 0.160243 0.366476 20° 0.16
58 0.4 0.65 0.16 1.0 0.349066 0.2 0.158335 0.354346 20° 0.16
59 0.4 0.65 0.16 1.0 0.349066 0.2 0.148304 0.235357 20° 0.16
60 0.4 0.65 0.16 1.0 0.349066 0.2 0.160054 0.206858 20° 0.16

Visualize 3D surfaces#

Hide definition: plot_mesh3d()

 
def plot_mesh3d(vertices, faces, *, title="", opacity=0.5, color="lightblue"):
    """Plot a triangular mesh with plotly."""
    I, J, K = faces.T
    x, y, z = vertices.T

    fig = go.Figure(
        data=[
            go.Mesh3d(
                x=x, y=y, z=z,
                i=I, j=J, k=K,
                opacity=opacity,
                color=color,
                showscale=False,
            ),
        ]
    )
    fig.update_layout(
        title=title,
        width=700, height=500,
        scene=dict(aspectmode="data"),
    )
    return fig

representative_ids = [0, 10, 20, 30, 40, 50]
bending_order = [str(deg) + "\u00b0" for deg in [20, 80, 140]]

for sid in representative_ids:
    label = (
        f"Specimen {sid + 1} "
        f"(bending={df_meta.iloc[sid]['bending_angle']}, "
        f"\u03b1={df_meta.iloc[sid]['aspect_ratio']})"
    )
    fig = plot_mesh3d(
        data.vertices[sid], data.faces[sid], title=label
    )
    fig.show()

Parameter mesh (disk parameterization)#

The parameter mesh maps each surface vertex to a point on the unit disk.

fig = plot_mesh3d(
    data.param_vertices[0], data.param_faces[0],
    title="Parameter mesh (specimen 1)", opacity=0.3,
)
fig.show()

Prepare disk parameterization#

DHA requires polar coordinates (r, theta) on the unit disk. We extract the (x, y) coordinates from the parameter mesh and convert them using xy2polar.

r_theta = [
    xy2polar(data.param_vertices[i][:, :2], centered=True)
    for i in range(len(data.vertices))
]

print(f"r_theta[0].shape: {r_theta[0].shape}")
print(f"r range: [{r_theta[0][:, 0].min():.4f}, {r_theta[0][:, 0].max():.4f}]")
print(f"theta range: [{r_theta[0][:, 1].min():.4f}, {r_theta[0][:, 1].max():.4f}]")

r_theta[0].shape: (8099, 2)
r range: [0.0066, 1.0000]
theta range: [0.0006, 6.2821]

X = data.vertices

Disk Harmonic Analysis#

dha = DiskHarmonicAnalysis(n_harmonics=15, n_dim=3)

coef = dha.fit_transform(X, r_theta=r_theta)
coef.shape

(60, 768)

Reconstruction from coefficients#

Reconstructing with different numbers of harmonics shows how DHA captures shape at different levels of detail. The inverse_transform method returns a regular (n_theta, n_r, 3) grid for each specimen.

sid = 0
I, J, K = data.faces[sid].T
x_orig, y_orig, z_orig = data.vertices[sid].T

n_max_values = [5, 10, 15]

for n_max in n_max_values:
    surf = dha.inverse_transform(coef[sid:sid+1], n_max=n_max)
    x_r, y_r, z_r = surf[0].transpose(2, 0, 1)

    fig = go.Figure(
        data=[
            go.Mesh3d(
                x=x_orig, y=y_orig, z=z_orig,
                i=I, j=J, k=K,
                opacity=0.3,
                color="lightblue",
                showscale=False,
                name="original",
            ),
            go.Surface(
                x=x_r, y=y_r, z=z_r,
                opacity=0.5,
                showscale=False,
                name=f"n_max={n_max}",
            ),
        ]
    )
    fig.update_layout(
        title=f"Reconstruction with n_max={n_max}",
        width=700, height=500,
        scene=dict(aspectmode="data"),
    )
    fig.show()