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.29720027e-02, 2.48897303e-03, 1.10721108e-02,
6.58749779e-02, -1.88412735e-03, -9.04224775e-04,
-2.13759717e-04, -3.48942457e-03, 1.14882134e-03,
8.96017559e-04, 7.86471071e-04, -4.42645170e-04,
-3.46744930e-04, 3.24924515e-04, 9.91263874e-04,
-1.96090736e-04, -7.51414278e-04, -2.44784490e-04,
4.44948628e-04, -1.12307100e-04, 2.26480956e-04,
5.49308705e-04, -5.92965002e-04, -3.68909956e-04,
-7.16184975e-04, 2.02576837e-04, -6.69970345e-04,
-1.32382675e-05, 1.73828508e-04, -9.36412973e-05,
-2.94812790e-04, 1.03155814e-03, 2.63623321e-04,
-5.80006547e-04, 1.01434629e-04, -8.53467267e-05,
3.10549552e-05, 1.14528700e-04, 5.37492303e-04,
-9.91802394e-06, 2.00438731e-04, 8.49202316e-05,
4.95629189e-05, 3.85075987e-04, 7.98220419e-06,
2.56847860e-04, 8.43471488e-06, -6.29500968e-04,
-3.00995302e-04, 7.15585807e-05, -1.84408001e-04,
-6.79017105e-05, 4.82502531e-05, 1.04724573e-06,
4.32939817e-05, 1.01194965e-04, 3.30893167e-06,
-1.10895371e-04, -5.50610205e-06, 2.64687831e-04,
-1.07519794e-04, -4.30962772e-04, -2.90266638e-05,
5.10077245e-05, -9.54234056e-05, -3.00744168e-05,
-1.92703580e-04, -1.39380619e-04, 1.70239138e-04,
1.55171285e-04, -2.37209124e-05, 4.18032759e-05,
-3.38266968e-05, -1.47348370e-04, -1.98609553e-05,
7.05419468e-05, 1.60770275e-05, -2.90684419e-05,
5.63951864e-05, 1.26833543e-04, -1.32156953e-06,
1.41823923e-04, 5.59828003e-05, 3.07813028e-06,
1.09042336e-04, -8.76470065e-06, -1.32338230e-05,
2.44472450e-05, -3.15469815e-05, -9.82344871e-08,
-9.65559342e-05, 5.18039336e-05, 4.19000719e-05,
-1.26956054e-04, 6.35386159e-05, 1.08406684e-04,
-8.52131926e-05, -1.00843301e-04, -1.27415102e-04,
-2.72490547e-05, -8.18067006e-05, -2.18692102e-05,
-1.80777047e-05, -7.16985716e-05, -7.77004296e-05,
4.37230935e-05, -5.84898377e-05, -1.86098273e-05,
-2.87965764e-05, -2.25188788e-05, -6.33736460e-05,
1.65206414e-04, 2.65670621e-05, -6.58327381e-05,
6.82958383e-05, 2.07534075e-05, -1.14730884e-04,
-1.78375467e-05, 2.24814396e-05, -1.72531366e-05,
-5.68174521e-05, 1.20166975e-01, -4.57018843e-02,
-2.92836043e-02, -1.78043164e-03, 2.08948481e-04,
1.17996547e-03, 4.58999547e-04, 6.39358734e-05,
1.08406660e-03, -2.43563367e-03, 6.71051521e-05,
-6.00420483e-04, 6.45703903e-05, 8.20336314e-05,
3.22800233e-04, 4.98165654e-04, -1.26058583e-03,
-1.02549119e-03, 1.29034852e-04, -1.08432119e-04,
8.58986048e-05, 2.80544222e-04, -5.41581392e-04,
-4.04062213e-04, -5.52273312e-04, 2.07193759e-04,
-5.98608185e-04, -7.48400072e-04, -2.41370407e-04,
-1.10436129e-04, -7.51653257e-07, 1.38446363e-04,
-4.02500974e-04, -3.91873261e-04, 3.25145353e-04,
2.61207831e-04, 2.90791749e-04, 6.13732496e-04,
4.02343839e-04, -3.18754783e-04, 1.88883125e-05,
-5.38619121e-04, 1.61558986e-04, -1.07760755e-04,
-2.59938628e-04, -2.28925063e-05, 3.48788585e-04,
8.19317590e-06, -9.31128809e-06, -1.06410469e-04,
-5.04272434e-05, 3.47659733e-04, -2.84962024e-05,
-1.59913824e-05, 1.76309126e-04, -2.37360878e-04,
4.02438848e-05, -1.45095609e-04, -9.18538735e-05,
-1.03949488e-04, 1.21149919e-04, 1.97081806e-04,
-1.18255381e-04, -6.77487977e-05, -7.69819396e-05,
5.34472158e-05, 2.21895992e-04, 3.02359660e-05,
-4.62163363e-05, 5.17774228e-05, -2.07755988e-05,
3.16721541e-05, 3.91490812e-05, -7.92975598e-05,
2.26871941e-04, -1.71638743e-05, 1.86382935e-04,
2.33975429e-04, 2.41957946e-04, -8.93480519e-05,
-8.02988059e-05, 1.01039714e-04, -1.35621396e-04,
1.13120542e-05, 7.80795546e-05, -7.84738943e-06,
-6.50121004e-05, 9.39321868e-05, 3.58259737e-06,
-3.36737006e-06, -1.16021175e-04, 1.09754137e-04,
5.83097921e-05, 7.76713433e-08, 2.95108811e-05,
1.66244394e-04, -1.59354817e-04, -2.87446623e-04,
1.87457041e-05, 2.14879360e-04, -4.45393524e-05,
-1.62349539e-06, -7.54175023e-05, -1.76685953e-04,
2.70132149e-05, -6.85878516e-05, -1.98476041e-05,
-2.33349579e-05, -9.86705287e-06, -4.99732090e-05,
-3.00361095e-05, -1.51402227e-05, -1.26094464e-05,
3.11205347e-05, -1.04432555e-05, 1.86333536e-05,
-1.26776861e-04, -9.79741865e-05, 4.12407581e-05,
-1.98545090e-06, -7.91947524e-05, -1.15047627e+00,
1.48890454e-02, -4.76763241e-02, 1.90055413e-03,
-3.25398189e-05, -1.71732166e-03, -2.83188382e-03,
3.97352567e-05, -9.23357260e-04, 1.16393902e-03,
2.99496985e-04, -2.07685095e-03, -1.95999273e-03,
-1.79953940e-03, -1.66321354e-03, 7.48336481e-05,
4.61072609e-04, 2.51663459e-04, 4.04768284e-04,
2.05540862e-05, -1.23425223e-03, 2.33905518e-04,
1.59711414e-04, -8.63769526e-05, 7.56225171e-04,
-1.16425726e-04, 3.57394244e-05, -8.97827418e-05,
2.08667656e-04, 2.49988843e-04, -2.23306819e-04,
4.02295913e-04, -2.86371128e-05, -2.39005403e-04,
3.29830528e-04, 8.69183654e-05, -2.99695953e-04,
-2.31110655e-04, -1.28784760e-04, -4.85688331e-04,
-7.46369608e-05, 7.29466595e-05, 4.07137608e-05,
-1.79749832e-04, -4.04248295e-04, 8.70209420e-05,
2.54030940e-04, 1.97306841e-04, -1.82018587e-04,
1.45122604e-04, -1.31216489e-05, -4.43416445e-05,
-1.98245894e-04, -7.64532009e-05, -1.27541228e-04,
1.71050351e-04, 1.05310585e-04, 5.84795483e-05,
-1.26664204e-04, 6.20618576e-05, 1.79014064e-04,
-1.06013529e-04, 1.29491545e-04, 5.15163343e-05,
1.67682138e-05, 3.95248770e-05, 1.85547144e-04,
-1.89267225e-05, -1.34466851e-05, -1.37260778e-04,
-5.84829378e-05, -2.40626995e-04, 3.29446467e-06,
-6.93285389e-05, 1.20621519e-04, -9.37505876e-05,
-1.87627882e-04, 3.80125445e-05, -7.76577074e-05,
5.22593975e-06, 4.15875718e-05, -6.35006514e-05,
-1.06864073e-04, -8.62110249e-05, 1.62823806e-04,
1.84039425e-04, -1.67959261e-05, 7.14654753e-06,
9.98702381e-06, -6.99501921e-05, -1.04087668e-04,
1.60687843e-05, 2.56547300e-04, -1.29621256e-04,
-2.21146565e-04, 2.63752359e-05, -2.88686628e-05,
1.06601977e-04, 6.10283153e-05, -7.49353051e-05,
3.81499296e-05, 2.47948450e-05, 5.01497808e-05,
4.89607884e-05, 8.87421086e-05, 8.50700623e-05,
-5.30953763e-05, 1.41718376e-04, 5.90517918e-05,
7.53406644e-05, -4.48072429e-05, -4.42725302e-05,
1.17121937e-04, 1.63044842e-05, -2.05474919e-05,
-8.24988011e-05, -6.62420418e-05, 3.30956938e-05,
9.80959775e-05, 6.17558677e-05, 4.57492397e-05]])
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()