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], [xyz2spherical(arr_para)])
X_transform
array([[-2.29720027e-02-2.61567645e-23j, 4.65806436e-02+1.75996971e-03j,
1.10721108e-02+6.43534645e-18j, -4.65806436e-02+1.75996971e-03j,
8.12339361e-04-1.33227922e-03j, -2.46739577e-03-6.39383470e-04j,
-2.13759717e-04+8.64495930e-18j, 2.46739577e-03-6.39383470e-04j,
8.12339361e-04+1.33227922e-03j, -1.38657089e-04+6.33580092e-04j,
7.00929407e-04+5.56119028e-04j, 2.29756328e-04-3.12997401e-04j,
-3.46744930e-04-3.53095454e-18j, -2.29756328e-04-3.12997401e-04j,
7.00929407e-04-5.56119028e-04j, 1.38657089e-04+6.33580092e-04j,
-5.06419253e-04-5.31330132e-04j, -2.60858732e-04-1.73088773e-04j,
-4.19289574e-04+3.14626192e-04j, 3.88419910e-04-7.94131122e-05j,
2.26480956e-04-7.65044475e-18j, -3.88419910e-04-7.94131122e-05j,
-4.19289574e-04-3.14626192e-04j, 2.60858732e-04-1.73088773e-04j,
-5.06419253e-04+5.31330132e-04j, -6.03492492e-05+1.43243455e-04j,
7.17251140e-05-4.73740574e-04j, -4.10126562e-04-9.36086869e-06j,
1.86409838e-04+1.22915317e-04j, 7.29421757e-04-6.62143963e-05j,
-2.94812790e-04-1.47976109e-21j, -7.29421757e-04-6.62143963e-05j,
1.86409838e-04-1.22915317e-04j, 4.10126562e-04-9.36086869e-06j,
7.17251140e-05+4.73740574e-04j, 6.03492492e-05+1.43243455e-04j,
-2.12835819e-04+2.19591694e-05j, -4.45124403e-04+8.09840203e-05j,
5.96424409e-06+3.80064453e-04j, 1.81618864e-04-7.01310198e-06j,
5.64427071e-06+1.41731586e-04j, 2.72289842e-04+6.00476717e-05j,
4.95629189e-05+3.68628739e-18j, -2.72289842e-04+6.00476717e-05j,
5.64427071e-06-1.41731586e-04j, -1.81618864e-04-7.01310198e-06j,
5.96424409e-06-3.80064453e-04j, 4.45124403e-04+8.09840203e-05j,
-2.12835819e-04-2.19591694e-05j, 3.60679079e-05+5.05995577e-05j,
-2.05249508e-05-1.30396148e-04j, -3.04736699e-04-4.80137600e-05j,
-7.60279752e-05+3.41180812e-05j, 1.87162560e-04+7.40514560e-07j,
-3.89340210e-06+3.06134680e-05j, -7.84148690e-05+7.15556462e-05j,
3.30893167e-06+1.40946282e-18j, 7.84148690e-05+7.15556462e-05j,
-3.89340210e-06-3.06134680e-05j, -1.87162560e-04+7.40514560e-07j,
-7.60279752e-05-3.41180812e-05j, 3.04736699e-04-4.80137600e-05j,
-2.05249508e-05+1.30396148e-04j, -3.60679079e-05+5.05995577e-05j,
-9.34490780e-07-6.74745372e-05j, 8.96848585e-05-2.12658241e-05j,
3.98774187e-05-1.36262008e-04j, -2.05544924e-05-9.85569806e-05j,
1.13681752e-05+1.20377249e-04j, 4.98806890e-05+1.09722668e-04j,
-1.40438162e-05-1.67732180e-05j, -1.04191032e-04+2.95593799e-05j,
-3.38266968e-05-1.84314369e-18j, 1.04191032e-04+2.95593799e-05j,
-1.40438162e-05+1.67732180e-05j, -4.98806890e-05+1.09722668e-04j,
1.13681752e-05-1.20377249e-04j, 2.05544924e-05-9.85569806e-05j,
3.98774187e-05+1.36262008e-04j, -8.96848585e-05-2.12658241e-05j,
-9.34490780e-07+6.74745372e-05j, -1.92679913e-05+1.00284658e-04j,
-9.00960825e-05+3.95858177e-05j, -7.13069823e-05+2.17656679e-06j,
-6.02548264e-05+7.71045749e-05j, 7.66551017e-05-6.19757926e-06j,
4.49285862e-05-9.35772598e-06j, -8.97714865e-05+1.72868127e-05j,
2.96278250e-05-2.23070845e-05j, 3.66309127e-05-6.94622720e-08j,
-9.65559342e-05-2.16840434e-19j, -3.66309127e-05-6.94622720e-08j,
2.96278250e-05+2.23070845e-05j, 8.97714865e-05+1.72868127e-05j,
4.49285862e-05+9.35772598e-06j, -7.66551017e-05-6.19757926e-06j,
-6.02548264e-05-7.71045749e-05j, 7.13069823e-05+2.17656679e-06j,
-9.00960825e-05-3.95858177e-05j, 1.92679913e-05+1.00284658e-04j,
-4.01760057e-05-5.78460728e-05j, -1.21998099e-05-1.54638669e-05j,
1.58967784e-05-1.27828676e-05j, -1.26130502e-05-5.06985462e-05j,
-8.11269864e-05-5.49425007e-05j, 1.46748752e-05+3.09168959e-05j,
4.82924504e-05-4.13585609e-05j, -4.65507756e-05-1.31591351e-05j,
1.87857498e-05-2.03622544e-05j, 1.16818575e-04-1.59232519e-05j,
-6.33736460e-05+1.30104261e-18j, -1.16818575e-04-1.59232519e-05j,
1.87857498e-05+2.03622544e-05j, 4.65507756e-05-1.31591351e-05j,
4.82924504e-05+4.13585609e-05j, -1.46748752e-05+3.09168959e-05j,
-8.11269864e-05+5.49425007e-05j, 1.26130502e-05-5.06985462e-05j,
1.58967784e-05+1.27828676e-05j, 1.21998099e-05-1.54638669e-05j,
-4.01760057e-05+5.78460728e-05j, 1.20166975e-01+3.89128987e-22j,
-1.25895529e-03-3.23161123e-02j, -2.92836043e-02+8.75152072e-19j,
1.25895529e-03-3.23161123e-02j, 7.66550844e-04+1.47748888e-04j,
4.52094896e-05+8.34361582e-04j, 4.58999547e-04-7.14178517e-18j,
-4.52094896e-05+8.34361582e-04j, 7.66550844e-04-1.47748888e-04j,
3.52256312e-04-1.72225309e-03j, 2.28254234e-04+4.74505081e-05j,
5.80065371e-05-4.24561395e-04j, 6.45703903e-05-2.28163855e-19j,
-5.80065371e-05-4.24561395e-04j, 2.28254234e-04-4.74505081e-05j,
-3.52256312e-04-1.72225309e-03j, -3.90516204e-04-8.91368787e-04j,
-2.85715131e-04-7.25131777e-04j, -3.82955875e-04+9.12414190e-05j,
1.98374721e-04-7.66730867e-05j, 8.58986048e-05-8.41942868e-18j,
-1.98374721e-04-7.66730867e-05j, -3.82955875e-04-9.12414190e-05j,
2.85715131e-04-7.25131777e-04j, -3.90516204e-04+8.91368787e-04j,
1.84701829e-04+1.46508112e-04j, 2.29912484e-04-4.23279907e-04j,
-2.77096240e-04-5.29198766e-04j, -2.84611168e-04-1.70674651e-04j,
9.78963623e-05-7.80901360e-05j, -7.51653257e-07-6.25107219e-18j,
-9.78963623e-05-7.80901360e-05j, -2.84611168e-04+1.70674651e-04j,
2.77096240e-04-5.29198766e-04j, 2.29912484e-04+4.23279907e-04j,
-1.84701829e-04+1.46508112e-04j, -6.58407495e-06+2.05620818e-04j,
5.79345024e-06+4.33974410e-04j, 2.46630773e-04+2.84500057e-04j,
-1.61874465e-05-2.25393668e-04j, -1.83804367e-04+1.33560538e-05j,
-7.61983608e-05-3.80861233e-04j, 1.61558986e-04-5.20417043e-18j,
7.61983608e-05-3.80861233e-04j, -1.83804367e-04-1.33560538e-05j,
1.61874465e-05-2.25393668e-04j, 2.46630773e-04-2.84500057e-04j,
-5.79345024e-06+4.33974410e-04j, -6.58407495e-06-2.05620818e-04j,
-4.79056343e-05-7.52435641e-05j, -8.36191816e-05-3.56574457e-05j,
1.39357882e-04+2.45832555e-04j, 8.56659296e-05-2.01498580e-05j,
-7.35033876e-05-1.13076150e-05j, -6.49504969e-05+1.24669379e-04j,
-1.02598089e-04-1.67839486e-04j, 4.02438848e-05-4.33680869e-19j,
1.02598089e-04-1.67839486e-04j, -6.49504969e-05-1.24669379e-04j,
7.35033876e-05-1.13076150e-05j, 8.56659296e-05+2.01498580e-05j,
-1.39357882e-04+2.45832555e-04j, -8.36191816e-05+3.56574457e-05j,
4.79056343e-05-7.52435641e-05j, -5.67798302e-05-5.44344515e-05j,
-6.31786134e-05+3.77928888e-05j, 1.71090104e-04+1.56904161e-04j,
1.65445612e-04+2.13800566e-05j, 1.31792637e-04-3.26798848e-05j,
-1.21366919e-05+3.66121668e-05j, 1.60422688e-04-1.46905668e-05j,
-5.60718422e-05+2.23955949e-05j, 3.91490812e-05-1.21430643e-17j,
5.60718422e-05+2.23955949e-05j, 1.60422688e-04+1.46905668e-05j,
1.21366919e-05+3.66121668e-05j, 1.31792637e-04+3.26798848e-05j,
-1.65445612e-04+2.13800566e-05j, 1.71090104e-04-1.56904161e-04j,
6.31786134e-05+3.77928888e-05j, -5.67798302e-05+5.44344515e-05j,
1.51942653e-04+7.14458667e-05j, 1.32552145e-05-9.58988085e-05j,
-2.03255457e-04+7.99883025e-06j, -1.12680872e-04+5.52105825e-05j,
1.17552538e-04-5.54894228e-06j, 2.08673441e-05-4.59704970e-05j,
5.49219335e-08+6.64200863e-05j, 4.12312494e-05+2.53327889e-06j,
7.76078944e-05-2.38109021e-06j, -1.16021175e-04-1.21430643e-17j,
-7.76078944e-05-2.38109021e-06j, 4.12312494e-05-2.53327889e-06j,
-5.49219336e-08+6.64200863e-05j, 2.08673441e-05+4.59704970e-05j,
-1.17552538e-04-5.54894228e-06j, -1.12680872e-04-5.52105825e-05j,
2.03255457e-04+7.99883025e-06j, 1.32552145e-05+9.58988085e-05j,
-1.51942653e-04+7.14458667e-05j, -5.59991464e-05-3.14940781e-05j,
-1.40392580e-06-1.14798460e-06j, 2.91616197e-05-5.33282273e-05j,
-6.92782117e-05-1.24935836e-04j, -8.96447780e-05+1.91012275e-05j,
1.31757707e-05-4.84989350e-05j, -7.38449680e-06-1.40343755e-05j,
2.20055411e-05-1.65003070e-05j, -8.91622508e-06-6.97705999e-06j,
-1.07057542e-05-3.53363950e-05j, -3.00361095e-05-5.20417043e-18j,
1.07057542e-05-3.53363950e-05j, -8.91622508e-06+6.97705999e-06j,
-2.20055411e-05-1.65003070e-05j, -7.38449680e-06+1.40343755e-05j,
-1.31757707e-05-4.84989350e-05j, -8.96447780e-05-1.91012275e-05j,
6.92782117e-05-1.24935836e-04j, 2.91616197e-05+5.33282273e-05j,
1.40392580e-06-1.14798460e-06j, -5.59991464e-05+3.14940781e-05j,
-1.15047627e+00-1.77800547e-22j, 1.34389471e-03+1.05281449e-02j,
-4.76763241e-02+6.00106604e-18j, -1.34389471e-03+1.05281449e-02j,
-6.52912180e-04-2.30091266e-05j, 2.80970695e-05-1.21432979e-03j,
-2.83188382e-03-3.22564718e-18j, -2.80970695e-05-1.21432979e-03j,
-6.52912180e-04+2.30091266e-05j, 5.29153800e-05+8.23029172e-04j,
-1.17606957e-03+2.11776349e-04j, -1.27246652e-03-1.46855539e-03j,
-1.95999273e-03+7.58274634e-18j, 1.27246652e-03-1.46855539e-03j,
-1.17606957e-03-2.11776349e-04j, -5.29153800e-05+8.23029172e-04j,
5.34731946e-04+3.26027569e-04j, -6.10777289e-05+1.77952939e-04j,
1.12933024e-04+2.86214398e-04j, 1.65396178e-04+1.45339338e-05j,
-1.23425223e-03+2.37316008e-18j, -1.65396178e-04+1.45339338e-05j,
1.12933024e-04-2.86214398e-04j, 6.10777289e-05+1.77952939e-04j,
5.34731946e-04-3.26027569e-04j, 6.14605656e-05-8.23254206e-05j,
2.33225403e-04+2.52715893e-05j, -1.69002341e-04-6.34859856e-05j,
-2.02494967e-05+1.47550315e-04j, 2.84466168e-04+1.76768806e-04j,
-2.23306819e-04+1.53150005e-18j, -2.84466168e-04+1.76768806e-04j,
-2.02494967e-05-1.47550315e-04j, 1.69002341e-04-6.34859856e-05j,
2.33225403e-04-2.52715893e-05j, -6.14605656e-05-8.23254206e-05j,
-1.28706577e-04-2.11917041e-04j, 1.39517005e-04-1.63419911e-04j,
1.79627001e-04-9.10645769e-05j, 6.15330982e-05-3.43433512e-04j,
-2.85846710e-04-5.27763011e-05j, -1.27102325e-04+5.15810776e-05j,
4.07137608e-05+8.67361738e-19j, 1.27102325e-04+5.15810776e-05j,
-2.85846710e-04+5.27763011e-05j, -6.15330982e-05-3.43433512e-04j,
1.79627001e-04+9.10645769e-05j, -1.39517005e-04-1.63419911e-04j,
-1.28706577e-04+2.11917041e-04j, 3.64275493e-05+1.02617177e-04j,
9.15643494e-05-9.27840694e-06j, -7.49628851e-05-3.13542775e-05j,
1.26582058e-04-1.40181016e-04j, 4.38843604e-05-5.40605768e-05j,
-8.95651174e-05-9.01852672e-05j, 4.13512852e-05+1.20950863e-04j,
1.05310585e-04+3.46944695e-18j, -4.13512852e-05+1.20950863e-04j,
-8.95651174e-05+9.01852672e-05j, -4.38843604e-05-5.40605768e-05j,
1.26582058e-04+1.40181016e-04j, 7.49628851e-05-3.13542775e-05j,
9.15643494e-05+9.27840694e-06j, -3.64275493e-05+1.02617177e-04j,
2.94068540e-05+1.18569177e-05j, 3.69529743e-06+2.79483085e-05j,
-5.49122915e-05+1.31201644e-04j, 2.68789280e-05-1.33832138e-05j,
-1.32672947e-04-9.50824223e-06j, -6.62916762e-05-9.70580269e-05j,
8.52922944e-05-4.13536819e-05j, -4.90226800e-05-1.70148980e-04j,
3.29446467e-06+1.73472348e-18j, 4.90226800e-05-1.70148980e-04j,
8.52922944e-05+4.13536819e-05j, 6.62916762e-05-9.70580269e-05j,
-1.32672947e-04+9.50824223e-06j, -2.68789280e-05-1.33832138e-05j,
-5.49122915e-05-1.31201644e-04j, -3.69529743e-06+2.79483085e-05j,
2.94068540e-05-1.18569177e-05j, -5.29872624e-05-4.49017412e-05j,
4.31535356e-05-7.55643108e-05j, 7.53789810e-05-6.09604003e-05j,
-2.04132272e-05+1.15133817e-04j, 1.86501082e-05+1.30135525e-04j,
-1.56374236e-04-1.18765133e-05j, -9.16560693e-05+5.05337222e-06j,
1.81406335e-04+7.06189226e-06j, 1.13623463e-05-4.94622552e-05j,
-1.04087668e-04+3.46944695e-18j, -1.13623463e-05-4.94622552e-05j,
1.81406335e-04-7.06189226e-06j, 9.16560693e-05+5.05337222e-06j,
-1.56374236e-04+1.18765133e-05j, -1.86501082e-05+1.30135525e-04j,
-2.04132272e-05-1.15133817e-04j, -7.53789810e-05-6.09604003e-05j,
4.31535356e-05+7.55643108e-05j, 5.29872624e-05-4.49017412e-05j,
3.23495976e-05+2.69760739e-05j, 4.36679928e-05+1.75326030e-05j,
6.93643309e-05+3.54612501e-05j, 2.34021895e-05+3.46205055e-05j,
-4.68401970e-05+6.27501468e-05j, -5.83354617e-05+6.01536180e-05j,
-1.45292709e-05-3.75441006e-05j, 1.15290114e-05+1.00210025e-04j,
8.28177162e-05+4.17559224e-05j, -3.13054063e-05+5.32738947e-05j,
-4.48072429e-05+2.16840434e-19j, 3.13054063e-05+5.32738947e-05j,
8.28177162e-05-4.17559224e-05j, -1.15290114e-05+1.00210025e-04j,
-1.45292709e-05+3.75441006e-05j, 5.83354617e-05+6.01536180e-05j,
-4.68401970e-05-6.27501468e-05j, -2.34021895e-05+3.46205055e-05j,
6.93643309e-05-3.54612501e-05j, -4.36679928e-05+1.75326030e-05j,
3.23495976e-05-2.69760739e-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()