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()