{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4834c5de",
   "metadata": {},
   "source": [
    "# Spherical Harmonic (SPHARM) Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35a45bc3",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "import shutil\n",
    "import urllib\n",
    "import tempfile\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "import plotly.graph_objects as go\n",
    "\n",
    "import vtk\n",
    "from vtk.numpy_interface import dataset_adapter as dsa\n",
    "\n",
    "from ktch.harmonic import SphericalHarmonicAnalysis, xyz2spherical"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a605bf4",
   "metadata": {},
   "source": [
    "## Load 3D potato surface data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6148550",
   "metadata": {},
   "outputs": [],
   "source": [
    "# parameter\n",
    "with urllib.request.urlopen(\n",
    "    \"https://strata.morphometrics.jp/examples/andesred_07_allSegments_para.vtk\"\n",
    ") as response:\n",
    "    with tempfile.NamedTemporaryFile() as tmp_file:\n",
    "        shutil.copyfileobj(response, tmp_file)\n",
    "        reader = vtk.vtkDataSetReader()\n",
    "        reader.SetFileName(tmp_file.name)\n",
    "        reader.Update()\n",
    "\n",
    "        dataset = reader.GetOutput()\n",
    "        obj_para = dsa.WrapDataObject(dataset)\n",
    "\n",
    "# surface\n",
    "with urllib.request.urlopen(\n",
    "    \"https://strata.morphometrics.jp/examples/andesred_07_allSegments_surf.vtk\"\n",
    ") as response:\n",
    "    with tempfile.NamedTemporaryFile() as tmp_file:\n",
    "        shutil.copyfileobj(response, tmp_file)\n",
    "        reader = vtk.vtkDataSetReader()\n",
    "        reader.SetFileName(tmp_file.name)\n",
    "        reader.Update()\n",
    "\n",
    "        dataset = reader.GetOutput()\n",
    "        obj_surf = dsa.WrapDataObject(dataset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "082c1773",
   "metadata": {},
   "outputs": [],
   "source": [
    "arr_para_faces = np.array(obj_para.Polygons.reshape(-1, 4)[:, 1:])\n",
    "arr_para = np.array(obj_para.Points)\n",
    "\n",
    "arr_surf_faces = np.array(obj_surf.Polygons.reshape(-1, 4)[:, 1:])\n",
    "arr_surf = np.array(obj_surf.Points)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94ef19bc",
   "metadata": {},
   "source": [
    "### Parameter mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f0fe821",
   "metadata": {},
   "outputs": [],
   "source": [
    "I, J, K = arr_para_faces.T\n",
    "x_surf, y_surf, z_surf = arr_para.T\n",
    "\n",
    "fig = go.Figure(\n",
    "    data=[\n",
    "        go.Mesh3d(\n",
    "            x=x_surf,\n",
    "            y=y_surf,\n",
    "            z=z_surf,\n",
    "            i=I,\n",
    "            j=J,\n",
    "            k=K,\n",
    "            opacity=0.3,\n",
    "            showscale=False,\n",
    "        ),\n",
    "    ]\n",
    ")\n",
    "\n",
    "fig.update_layout(\n",
    "    width=700,\n",
    "    height=700,\n",
    "    autosize=False,\n",
    "    scene=dict(\n",
    "        camera=dict(\n",
    "            up=dict(x=0, y=0, z=1),\n",
    "            eye=dict(\n",
    "                x=1.1,\n",
    "                y=1.1,\n",
    "                z=1.1,\n",
    "            ),\n",
    "        ),\n",
    "        aspectmode=\"data\",\n",
    "    ),\n",
    ")\n",
    "\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95571a7f",
   "metadata": {},
   "source": [
    "### Surface mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68573eff",
   "metadata": {},
   "outputs": [],
   "source": [
    "I, J, K = arr_surf_faces.T\n",
    "x_surf, y_surf, z_surf = arr_surf.T\n",
    "\n",
    "fig = go.Figure(\n",
    "    data=[\n",
    "        go.Mesh3d(\n",
    "            x=x_surf,\n",
    "            y=y_surf,\n",
    "            z=z_surf,\n",
    "            i=I,\n",
    "            j=J,\n",
    "            k=K,\n",
    "            opacity=0.3,\n",
    "            showscale=False,\n",
    "        ),\n",
    "    ]\n",
    ")\n",
    "\n",
    "fig.update_layout(\n",
    "    width=700,\n",
    "    height=700,\n",
    "    autosize=False,\n",
    "    scene=dict(\n",
    "        camera=dict(\n",
    "            up=dict(x=0, y=0, z=1),\n",
    "            eye=dict(\n",
    "                x=1.1,\n",
    "                y=1.1,\n",
    "                z=1.1,\n",
    "            ),\n",
    "        ),\n",
    "        aspectmode=\"data\",\n",
    "    ),\n",
    ")\n",
    "\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c0fd13e",
   "metadata": {},
   "source": [
    "## SPHARM Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a85f2ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "spharm_analysis = SphericalHarmonicAnalysis(n_harmonics=10)\n",
    "X_transform = spharm_analysis.fit_transform([arr_surf], [xyz2spherical(arr_para)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd57b9fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_transform"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d5a3a20",
   "metadata": {},
   "source": [
    "## Inverse analysis\n",
    "\n",
    "Reconstruct 3D surface shapes from SPHARM coefficients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24921b8a",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_coords = spharm_analysis.inverse_transform(X_transform.reshape(1, 3, -1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9583c671",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_sph, y_sph, z_sph = np.real(X_coords[0].T)\n",
    "\n",
    "fig = go.Figure(\n",
    "    data=[\n",
    "        go.Surface(\n",
    "            x=x_sph,\n",
    "            y=y_sph,\n",
    "            z=z_sph,\n",
    "            opacity=0.3,\n",
    "            showscale=False,\n",
    "        ),\n",
    "    ]\n",
    ")\n",
    "\n",
    "fig.update_layout(\n",
    "    width=700,\n",
    "    height=700,\n",
    "    autosize=False,\n",
    "    scene=dict(\n",
    "        camera=dict(\n",
    "            up=dict(x=0, y=0, z=1),\n",
    "            eye=dict(\n",
    "                x=1.1,\n",
    "                y=1.1,\n",
    "                z=1.1,\n",
    "            ),\n",
    "        ),\n",
    "        aspectmode=\"data\",\n",
    "    ),\n",
    ")\n",
    "\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0367c0bf",
   "metadata": {},
   "source": [
    "### Original vs reconstructed shapes\n",
    "Overlay the original and reconstructed surface data.\n",
    "This demonstrates that the original shape is well reproduced."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5c228ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "I, J, K = arr_surf_faces.T\n",
    "x_surf, y_surf, z_surf = arr_surf.T\n",
    "\n",
    "x_sph, y_sph, z_sph = np.real(X_coords[0].T)\n",
    "\n",
    "\n",
    "fig = go.Figure(\n",
    "    data=[\n",
    "        go.Mesh3d(\n",
    "            x=x_surf,\n",
    "            y=y_surf,\n",
    "            z=z_surf,\n",
    "            i=I,\n",
    "            j=J,\n",
    "            k=K,\n",
    "            opacity=0.8,\n",
    "            showscale=False,\n",
    "        ),\n",
    "        go.Surface(\n",
    "            x=x_sph,\n",
    "            y=y_sph,\n",
    "            z=z_sph,\n",
    "            opacity=0.4,\n",
    "            showscale=False,\n",
    "        ),\n",
    "    ]\n",
    ")\n",
    "\n",
    "fig.update_layout(\n",
    "    width=700,\n",
    "    height=700,\n",
    "    autosize=False,\n",
    "    scene=dict(\n",
    "        camera=dict(\n",
    "            up=dict(x=0, y=0, z=1),\n",
    "            eye=dict(\n",
    "                x=1.1,\n",
    "                y=1.1,\n",
    "                z=1.1,\n",
    "            ),\n",
    "        ),\n",
    "        aspectmode=\"data\",\n",
    "    ),\n",
    ")\n",
    "\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a558b40",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "jupytext": {
   "default_lexer": "ipython3"
  },
  "kernelspec": {
   "display_name": "ktch",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}