Spherical Harmonic (SPHARM) Analysis#
from pathlib import Path
import shutil
import urllib.request
import tempfile
import numpy as np
import plotly.graph_objects as go
import vtk
from vtk.numpy_interface import dataset_adapter as dsa
from ktch.harmonic import SphericalHarmonicAnalysis, xyz2spherical
Load 3D potato surface data#
# parameter
req = urllib.request.Request(
"https://strata.morphometrics.jp/examples/danshaku_08_allSegments_para.vtp",
headers={"User-Agent": "Mozilla/5.0"},
)
with urllib.request.urlopen(req) as response:
with tempfile.NamedTemporaryFile() as tmp_file:
shutil.copyfileobj(response, tmp_file)
reader = vtk.vtkXMLPolyDataReader()
reader.SetFileName(tmp_file.name)
reader.Update()
dataset = reader.GetOutput()
obj_para = dsa.WrapDataObject(dataset)
# surface
req = urllib.request.Request(
"https://strata.morphometrics.jp/examples/danshaku_08_allSegments_surf.vtp",
headers={"User-Agent": "Mozilla/5.0"},
)
with urllib.request.urlopen(req) as response:
with tempfile.NamedTemporaryFile() as tmp_file:
shutil.copyfileobj(response, tmp_file)
reader = vtk.vtkXMLPolyDataReader()
reader.SetFileName(tmp_file.name)
reader.Update()
dataset = reader.GetOutput()
obj_surf = dsa.WrapDataObject(dataset)
arr_para_faces = np.array(obj_para.Polygons.reshape(-1, 4)[:, 1:])
arr_para = np.array(obj_para.Points)
arr_surf_faces = np.array(obj_surf.Polygons.reshape(-1, 4)[:, 1:])
arr_surf = np.array(obj_surf.Points)
Parameter mesh#
I, J, K = arr_para_faces.T
x_surf, y_surf, z_surf = arr_para.T
fig = go.Figure(
data=[
go.Mesh3d(
x=x_surf,
y=y_surf,
z=z_surf,
i=I,
j=J,
k=K,
opacity=0.3,
showscale=False,
),
]
)
fig.update_layout(
width=700,
height=700,
autosize=False,
scene=dict(
camera=dict(
up=dict(x=0, y=0, z=1),
eye=dict(
x=1.1,
y=1.1,
z=1.1,
),
),
aspectmode="data",
),
)
fig.show()
Surface mesh#
I, J, K = arr_surf_faces.T
x_surf, y_surf, z_surf = arr_surf.T
fig = go.Figure(
data=[
go.Mesh3d(
x=x_surf,
y=y_surf,
z=z_surf,
i=I,
j=J,
k=K,
opacity=0.3,
showscale=False,
),
]
)
fig.update_layout(
width=700,
height=700,
autosize=False,
scene=dict(
camera=dict(
up=dict(x=0, y=0, z=1),
eye=dict(
x=1.1,
y=1.1,
z=1.1,
),
),
aspectmode="data",
),
)
fig.show()
SPHARM Analysis#
spharm_analysis = SphericalHarmonicAnalysis(n_harmonics=10)
X_transform = spharm_analysis.fit_transform([arr_surf], theta_phi=[xyz2spherical(arr_para)])
X_transform
array([[-2.29719996e-02-3.97241051e-23j, 4.65806432e-02+1.75996968e-03j,
1.10721103e-02+2.84109639e-18j, -4.65806432e-02+1.75996968e-03j,
8.12338600e-04-1.33227925e-03j, -2.46739700e-03-6.39383547e-04j,
-2.13760644e-04+1.54404665e-17j, 2.46739700e-03-6.39383547e-04j,
8.12338600e-04+1.33227925e-03j, -1.38656782e-04+6.33580174e-04j,
7.00929321e-04+5.56119178e-04j, 2.29757474e-04-3.12997289e-04j,
-3.46744637e-04+1.84782594e-17j, -2.29757474e-04-3.12997289e-04j,
7.00929321e-04-5.56119178e-04j, 1.38656782e-04+6.33580174e-04j,
-5.06419329e-04-5.31330313e-04j, -2.60858784e-04-1.73088926e-04j,
-4.19289743e-04+3.14626066e-04j, 3.88419461e-04-7.94131477e-05j,
2.26480784e-04+8.03934580e-18j, -3.88419461e-04-7.94131477e-05j,
-4.19289743e-04-3.14626066e-04j, 2.60858784e-04-1.73088926e-04j,
-5.06419329e-04+5.31330313e-04j, -6.03488876e-05+1.43243633e-04j,
7.17251386e-05-4.73740456e-04j, -4.10126379e-04-9.36077164e-06j,
1.86409927e-04+1.22915339e-04j, 7.29422056e-04-6.62143981e-05j,
-2.94812660e-04-1.02641327e-17j, -7.29422056e-04-6.62143981e-05j,
1.86409927e-04-1.22915339e-04j, 4.10126379e-04-9.36077164e-06j,
7.17251386e-05+4.73740456e-04j, 6.03488876e-05+1.43243633e-04j,
-2.12835911e-04+2.19590277e-05j, -4.45124512e-04+8.09840314e-05j,
5.96413558e-06+3.80064498e-04j, 1.81618901e-04-7.01313295e-06j,
5.64420824e-06+1.41731599e-04j, 2.72289488e-04+6.00476602e-05j,
4.95628895e-05-3.25260652e-19j, -2.72289488e-04+6.00476602e-05j,
5.64420824e-06-1.41731599e-04j, -1.81618901e-04-7.01313295e-06j,
5.96413558e-06-3.80064498e-04j, 4.45124512e-04+8.09840314e-05j,
-2.12835911e-04-2.19590277e-05j, 3.60677760e-05+5.05996248e-05j,
-2.05248229e-05-1.30396228e-04j, -3.04736528e-04-4.80138500e-05j,
-7.60279686e-05+3.41180263e-05j, 1.87162552e-04+7.40527552e-07j,
-3.89338380e-06+3.06134786e-05j, -7.84148016e-05+7.15556321e-05j,
3.30879830e-06-5.42101086e-18j, 7.84148016e-05+7.15556321e-05j,
-3.89338380e-06-3.06134786e-05j, -1.87162552e-04+7.40527552e-07j,
-7.60279686e-05-3.41180263e-05j, 3.04736528e-04-4.80138500e-05j,
-2.05248229e-05+1.30396228e-04j, -3.60677760e-05+5.05996248e-05j,
-9.34574869e-07-6.74746145e-05j, 8.96847979e-05-2.12657599e-05j,
3.98774096e-05-1.36261971e-04j, -2.05544733e-05-9.85569776e-05j,
1.13681941e-05+1.20377279e-04j, 4.98806795e-05+1.09722649e-04j,
-1.40439009e-05-1.67731999e-05j, -1.04191151e-04+2.95593873e-05j,
-3.38266988e-05+5.20417043e-18j, 1.04191151e-04+2.95593873e-05j,
-1.40439009e-05+1.67731999e-05j, -4.98806795e-05+1.09722649e-04j,
1.13681941e-05-1.20377279e-04j, 2.05544733e-05-9.85569776e-05j,
3.98774096e-05+1.36261971e-04j, -8.96847979e-05-2.12657599e-05j,
-9.34574869e-07+6.74746145e-05j, -1.92679583e-05+1.00284763e-04j,
-9.00960878e-05+3.95857807e-05j, -7.13071262e-05+2.17659012e-06j,
-6.02548200e-05+7.71045946e-05j, 7.66551449e-05-6.19757743e-06j,
4.49285893e-05-9.35775164e-06j, -8.97714614e-05+1.72867830e-05j,
2.96278453e-05-2.23070864e-05j, 3.66311942e-05-6.94437548e-08j,
-9.65559405e-05+4.11996826e-18j, -3.66311942e-05-6.94437548e-08j,
2.96278453e-05+2.23070864e-05j, 8.97714614e-05+1.72867830e-05j,
4.49285893e-05+9.35775164e-06j, -7.66551449e-05-6.19757743e-06j,
-6.02548200e-05-7.71045946e-05j, 7.13071262e-05+2.17659012e-06j,
-9.00960878e-05-3.95857807e-05j, 1.92679583e-05+1.00284763e-04j,
-4.01760016e-05-5.78461739e-05j, -1.21997871e-05-1.54638554e-05j,
1.58968167e-05-1.27828614e-05j, -1.26130155e-05-5.06985809e-05j,
-8.11270015e-05-5.49425179e-05j, 1.46748589e-05+3.09168796e-05j,
4.82924408e-05-4.13585317e-05j, -4.65507421e-05-1.31591393e-05j,
1.87857438e-05-2.03622866e-05j, 1.16818448e-04-1.59232602e-05j,
-6.33736662e-05+2.16840434e-19j, -1.16818448e-04-1.59232602e-05j,
1.87857438e-05+2.03622866e-05j, 4.65507421e-05-1.31591393e-05j,
4.82924408e-05+4.13585317e-05j, -1.46748589e-05+3.09168796e-05j,
-8.11270015e-05+5.49425179e-05j, 1.26130155e-05-5.06985809e-05j,
1.58968167e-05+1.27828614e-05j, 1.21997871e-05-1.54638554e-05j,
-4.01760016e-05+5.78461739e-05j, 1.20166973e-01-1.64814179e-22j,
-1.25895529e-03-3.23161133e-02j, -2.92836034e-02+1.29811592e-17j,
1.25895529e-03-3.23161133e-02j, 7.66550960e-04+1.47750008e-04j,
4.52095771e-05+8.34362517e-04j, 4.58999997e-04-2.88512889e-18j,
-4.52095771e-05+8.34362517e-04j, 7.66550960e-04-1.47750008e-04j,
3.52256193e-04-1.72225418e-03j, 2.28254141e-04+4.74504799e-05j,
5.80065413e-05-4.24562173e-04j, 6.45702938e-05-2.31889492e-18j,
-5.80065413e-05-4.24562173e-04j, 2.28254141e-04-4.74504799e-05j,
-3.52256193e-04-1.72225418e-03j, -3.90516140e-04-8.91368221e-04j,
-2.85715181e-04-7.25131805e-04j, -3.82955943e-04+9.12415337e-05j,
1.98374714e-04-7.66727158e-05j, 8.58986943e-05+1.91122843e-18j,
-1.98374714e-04-7.66727158e-05j, -3.82955943e-04-9.12415337e-05j,
2.85715181e-04-7.25131805e-04j, -3.90516140e-04+8.91368221e-04j,
1.84701860e-04+1.46507852e-04j, 2.29912598e-04-4.23279707e-04j,
-2.77096138e-04-5.29199007e-04j, -2.84611149e-04-1.70674693e-04j,
9.78962874e-05-7.80904003e-05j, -7.51987747e-07+7.70480144e-18j,
-9.78962874e-05-7.80904003e-05j, -2.84611149e-04+1.70674693e-04j,
2.77096138e-04-5.29199007e-04j, 2.29912598e-04+4.23279707e-04j,
-1.84701860e-04+1.46507852e-04j, -6.58410738e-06+2.05621055e-04j,
5.79339377e-06+4.33974443e-04j, 2.46630728e-04+2.84500311e-04j,
-1.61874376e-05-2.25393588e-04j, -1.83804358e-04+1.33561875e-05j,
-7.61983544e-05-3.80861088e-04j, 1.61558962e-04-3.46944695e-18j,
7.61983544e-05-3.80861088e-04j, -1.83804358e-04-1.33561875e-05j,
1.61874376e-05-2.25393588e-04j, 2.46630728e-04-2.84500311e-04j,
-5.79339377e-06+4.33974443e-04j, -6.58410738e-06-2.05621055e-04j,
-4.79056681e-05-7.52437781e-05j, -8.36191517e-05-3.56575192e-05j,
1.39357933e-04+2.45832330e-04j, 8.56659637e-05-2.01499157e-05j,
-7.35033977e-05-1.13077138e-05j, -6.49504589e-05+1.24669381e-04j,
-1.02598071e-04-1.67839479e-04j, 4.02440729e-05-1.73472348e-18j,
1.02598071e-04-1.67839479e-04j, -6.49504589e-05-1.24669381e-04j,
7.35033977e-05-1.13077138e-05j, 8.56659637e-05+2.01499157e-05j,
-1.39357933e-04+2.45832330e-04j, -8.36191517e-05+3.56575192e-05j,
4.79056681e-05-7.52437781e-05j, -5.67797902e-05-5.44342159e-05j,
-6.31786434e-05+3.77929190e-05j, 1.71090071e-04+1.56904319e-04j,
1.65445586e-04+2.13800911e-05j, 1.31792650e-04-3.26798107e-05j,
-1.21366624e-05+3.66121517e-05j, 1.60422661e-04-1.46905327e-05j,
-5.60717948e-05+2.23956988e-05j, 3.91490777e-05+0.00000000e+00j,
5.60717948e-05+2.23956988e-05j, 1.60422661e-04+1.46905327e-05j,
1.21366624e-05+3.66121517e-05j, 1.31792650e-04+3.26798107e-05j,
-1.65445586e-04+2.13800911e-05j, 1.71090071e-04-1.56904319e-04j,
6.31786434e-05+3.77929190e-05j, -5.67797902e-05+5.44342159e-05j,
1.51942633e-04+7.14456428e-05j, 1.32552332e-05-9.58988705e-05j,
-2.03255482e-04+7.99868710e-06j, -1.12680959e-04+5.52105218e-05j,
1.17552493e-04-5.54903622e-06j, 2.08673148e-05-4.59704986e-05j,
5.49221255e-08+6.64200042e-05j, 4.12312102e-05+2.53328437e-06j,
7.76078590e-05-2.38128819e-06j, -1.16021144e-04-3.46944695e-18j,
-7.76078590e-05-2.38128819e-06j, 4.12312102e-05-2.53328437e-06j,
-5.49221255e-08+6.64200042e-05j, 2.08673148e-05+4.59704986e-05j,
-1.17552493e-04-5.54903622e-06j, -1.12680959e-04-5.52105218e-05j,
2.03255482e-04+7.99868710e-06j, 1.32552332e-05+9.58988705e-05j,
-1.51942633e-04+7.14456428e-05j, -5.59991485e-05-3.14939260e-05j,
-1.40394293e-06-1.14789166e-06j, 2.91616937e-05-5.33280856e-05j,
-6.92781531e-05-1.24935864e-04j, -8.96447539e-05+1.91012666e-05j,
1.31757350e-05-4.84989437e-05j, -7.38447758e-06-1.40342959e-05j,
2.20055472e-05-1.65003067e-05j, -8.91618487e-06-6.97705342e-06j,
-1.07057714e-05-3.53363045e-05j, -3.00361050e-05-6.07153217e-18j,
1.07057714e-05-3.53363045e-05j, -8.91618487e-06+6.97705342e-06j,
-2.20055472e-05-1.65003067e-05j, -7.38447758e-06+1.40342959e-05j,
-1.31757350e-05-4.84989437e-05j, -8.96447539e-05-1.91012666e-05j,
6.92781531e-05-1.24935864e-04j, 2.91616937e-05+5.33280856e-05j,
1.40394293e-06-1.14789166e-06j, -5.59991485e-05+3.14939260e-05j,
-1.15047627e+00+6.54215011e-22j, 1.34389471e-03+1.05281453e-02j,
-4.76763228e-02+4.83468792e-18j, -1.34389471e-03+1.05281453e-02j,
-6.52912285e-04-2.30094849e-05j, 2.80968562e-05-1.21433026e-03j,
-2.83188315e-03-2.25213608e-18j, -2.80968562e-05-1.21433026e-03j,
-6.52912285e-04+2.30094849e-05j, 5.29154744e-05+8.23029550e-04j,
-1.17606955e-03+2.11776435e-04j, -1.27246643e-03-1.46855510e-03j,
-1.95999305e-03-9.66776640e-18j, 1.27246643e-03-1.46855510e-03j,
-1.17606955e-03-2.11776435e-04j, -5.29154744e-05+8.23029550e-04j,
5.34731954e-04+3.26027356e-04j, -6.10776699e-05+1.77952898e-04j,
1.12933138e-04+2.86214395e-04j, 1.65396327e-04+1.45339725e-05j,
-1.23425201e-03+6.28357957e-18j, -1.65396327e-04+1.45339725e-05j,
1.12933138e-04-2.86214395e-04j, 6.10776699e-05+1.77952898e-04j,
5.34731954e-04-3.26027356e-04j, 6.14605455e-05-8.23253543e-05j,
2.33225376e-04+2.52715718e-05j, -1.69002335e-04-6.34859981e-05j,
-2.02495901e-05+1.47550310e-04j, 2.84466018e-04+1.76768787e-04j,
-2.23307368e-04-3.18155981e-18j, -2.84466018e-04+1.76768787e-04j,
-2.02495901e-05-1.47550310e-04j, 1.69002335e-04-6.34859981e-05j,
2.33225376e-04-2.52715718e-05j, -6.14605455e-05-8.23253543e-05j,
-1.28706568e-04-2.11917118e-04j, 1.39517002e-04-1.63419949e-04j,
1.79627045e-04-9.10646293e-05j, 6.15331434e-05-3.43433553e-04j,
-2.85846701e-04-5.27763340e-05j, -1.27102298e-04+5.15810600e-05j,
4.07137668e-05-4.55364912e-18j, 1.27102298e-04+5.15810600e-05j,
-2.85846701e-04+5.27763340e-05j, -6.15331434e-05-3.43433553e-04j,
1.79627045e-04+9.10646293e-05j, -1.39517002e-04-1.63419949e-04j,
-1.28706568e-04+2.11917118e-04j, 3.64275765e-05+1.02617308e-04j,
9.15643396e-05-9.27832127e-06j, -7.49629484e-05-3.13542099e-05j,
1.26582016e-04-1.40180969e-04j, 4.38843227e-05-5.40605431e-05j,
-8.95650925e-05-9.01852776e-05j, 4.13512733e-05+1.20950801e-04j,
1.05310824e-04-2.16840434e-18j, -4.13512733e-05+1.20950801e-04j,
-8.95650925e-05+9.01852776e-05j, -4.38843227e-05-5.40605431e-05j,
1.26582016e-04+1.40180969e-04j, 7.49629484e-05-3.13542099e-05j,
9.15643396e-05+9.27832127e-06j, -3.64275765e-05+1.02617308e-04j,
2.94068236e-05+1.18567677e-05j, 3.69529882e-06+2.79482451e-05j,
-5.49122756e-05+1.31201617e-04j, 2.68789428e-05-1.33831890e-05j,
-1.32672991e-04-9.50820777e-06j, -6.62916982e-05-9.70580634e-05j,
8.52923100e-05-4.13536739e-05j, -4.90226472e-05-1.70149051e-04j,
3.29447035e-06-5.58364119e-18j, 4.90226472e-05-1.70149051e-04j,
8.52923100e-05+4.13536739e-05j, 6.62916982e-05-9.70580634e-05j,
-1.32672991e-04+9.50820777e-06j, -2.68789428e-05-1.33831890e-05j,
-5.49122756e-05-1.31201617e-04j, -3.69529882e-06+2.79482451e-05j,
2.94068236e-05-1.18567677e-05j, -5.29872355e-05-4.49016416e-05j,
4.31535752e-05-7.55642781e-05j, 7.53789970e-05-6.09603945e-05j,
-2.04132347e-05+1.15133824e-04j, 1.86500814e-05+1.30135546e-04j,
-1.56374211e-04-1.18765169e-05j, -9.16560292e-05+5.05342137e-06j,
1.81406339e-04+7.06191673e-06j, 1.13623673e-05-4.94621053e-05j,
-1.04087596e-04+4.33680869e-18j, -1.13623673e-05-4.94621053e-05j,
1.81406339e-04-7.06191673e-06j, 9.16560292e-05+5.05342137e-06j,
-1.56374211e-04+1.18765169e-05j, -1.86500814e-05+1.30135546e-04j,
-2.04132347e-05-1.15133824e-04j, -7.53789970e-05-6.09603945e-05j,
4.31535752e-05+7.55642781e-05j, 5.29872355e-05-4.49016416e-05j,
3.23495888e-05+2.69759906e-05j, 4.36679765e-05+1.75325729e-05j,
6.93643026e-05+3.54612592e-05j, 2.34021982e-05+3.46204663e-05j,
-4.68401949e-05+6.27501050e-05j, -5.83354632e-05+6.01535955e-05j,
-1.45292659e-05-3.75441073e-05j, 1.15290101e-05+1.00209995e-04j,
8.28176897e-05+4.17559136e-05j, -3.13054105e-05+5.32738639e-05j,
-4.48072337e-05+2.60208521e-18j, 3.13054105e-05+5.32738639e-05j,
8.28176897e-05-4.17559136e-05j, -1.15290101e-05+1.00209995e-04j,
-1.45292659e-05+3.75441073e-05j, 5.83354632e-05+6.01535955e-05j,
-4.68401949e-05-6.27501050e-05j, -2.34021982e-05+3.46204663e-05j,
6.93643026e-05-3.54612592e-05j, -4.36679765e-05+1.75325729e-05j,
3.23495888e-05-2.69759906e-05j]])
Inverse analysis#
Reconstruct 3D surface shapes from SPHARM coefficients.
X_coords = spharm_analysis.inverse_transform(X_transform.reshape(1, 3, -1))
x_sph, y_sph, z_sph = np.real(X_coords[0].T)
fig = go.Figure(
data=[
go.Surface(
x=x_sph,
y=y_sph,
z=z_sph,
opacity=0.3,
showscale=False,
),
]
)
fig.update_layout(
width=700,
height=700,
autosize=False,
scene=dict(
camera=dict(
up=dict(x=0, y=0, z=1),
eye=dict(
x=1.1,
y=1.1,
z=1.1,
),
),
aspectmode="data",
),
)
fig.show()
Original vs reconstructed shapes#
Overlay the original and reconstructed surface data. This demonstrates that the original shape is well reproduced.
I, J, K = arr_surf_faces.T
x_surf, y_surf, z_surf = arr_surf.T
x_sph, y_sph, z_sph = np.real(X_coords[0].T)
fig = go.Figure(
data=[
go.Mesh3d(
x=x_surf,
y=y_surf,
z=z_surf,
i=I,
j=J,
k=K,
opacity=0.8,
showscale=False,
),
go.Surface(
x=x_sph,
y=y_sph,
z=z_sph,
opacity=0.4,
showscale=False,
),
]
)
fig.update_layout(
width=700,
height=700,
autosize=False,
scene=dict(
camera=dict(
up=dict(x=0, y=0, z=1),
eye=dict(
x=1.1,
y=1.1,
z=1.1,
),
),
aspectmode="data",
),
)
fig.show()