From 8fac11480cdd6b18fb9006787b4521ce3e2a99b1 Mon Sep 17 00:00:00 2001
From: Mariska Wesseling <M.G.H.Wesseling@tudelft.nl>
Date: Fri, 17 Jun 2022 17:16:01 +0200
Subject: [PATCH] updates to notebooks

---
 LigamentStudy/VisualiseSSM.py        |  47 +--
 LigamentStudy/VisualizeMeanSSM.ipynb | 183 +++++++++--
 LigamentStudy/testVisualizeSSM.py    | 441 +++++++++++++++++++++++++++
 3 files changed, 631 insertions(+), 40 deletions(-)
 create mode 100644 LigamentStudy/testVisualizeSSM.py

diff --git a/LigamentStudy/VisualiseSSM.py b/LigamentStudy/VisualiseSSM.py
index 2414586..6d64647 100644
--- a/LigamentStudy/VisualiseSSM.py
+++ b/LigamentStudy/VisualiseSSM.py
@@ -296,6 +296,7 @@ if __name__ == '__main__':
             mapper2.SetInputData(point_cloud_lig)
             actor2 = vtk.vtkActor()
             actor2.SetMapper(mapper2)
+            actor2.GetProperty().RenderPointsAsSpheresOn()
             actor2.GetProperty().SetColor(1, 0, 0)
             actor2.GetProperty().SetPointSize(7.5)
 
@@ -314,6 +315,7 @@ if __name__ == '__main__':
             surf_actor.GetProperty().SetOpacity(1.0)
 
             legend = vtk.vtkScalarBarActor()
+            legend.SetOrientationToHorizontal()
             labelFormat = vtk.vtkTextProperty()
             labelFormat.SetFontSize(16)
             titleFormat = vtk.vtkTextProperty()
@@ -338,10 +340,13 @@ if __name__ == '__main__':
             # text_prop_cb.SetFontSize(500)
             text_prop_cb.ShadowOff()
             legend.SetLabelTextProperty(text_prop_cb)
-            legend.SetMaximumWidthInPixels(75)
-            legend.SetMaximumHeightInPixels(300)
+            # legend.SetMaximumWidthInPixels(75)
+            # legend.SetMaximumHeightInPixels(300)
+            legend.SetMaximumWidthInPixels(300)
+            legend.SetMaximumHeightInPixels(75)
             legend.SetTitleTextProperty(text_prop_cb)
-            legend.SetPosition(0.85,0.5)
+            # legend.SetPosition(0.85,0.5)
+            legend.SetPosition(0.5, 0.85)
 
             # Renderer
             renderer = vtk.vtkRenderer()
@@ -349,7 +354,7 @@ if __name__ == '__main__':
             renderer.AddActor(actor2)
             renderer.AddActor(bone_actor)
             # renderer.AddActor(spline_actor)
-            renderer.AddActor(surf_actor)
+            # renderer.AddActor(surf_actor)
             if not subject == 100 and not subject == 'S0':
                 renderer.AddActor(wire_actor)
                 # renderer.AddActor(wire_actor2)
@@ -377,23 +382,23 @@ if __name__ == '__main__':
             renderWindowInteractor.Start()
 
 
-# polyData = vtk.vtkPolyData()
-# polyData.DeepCopy(surf_actor.GetMapper().GetInput())
-# transform = vtk.vtkTransform()
-# transform.SetMatrix(surf_actor.GetMatrix())
-# fil = vtk.vtkTransformPolyDataFilter()
-# fil.SetTransform(transform)
-# fil.SetInputDataObject(polyData)
-# fil.Update()
-# polyData.DeepCopy(fil.GetOutput())
-#
-# writer = vtk.vtkPLYWriter()
-# writer.SetFileTypeToASCII()
-# writer.SetColorModeToDefault()
-# filename = r'C:\Users\mariskawesseli\Documents\GitLab\tibia_surf_ply.ply'
-# writer.SetFileName(filename)
-# writer.SetInputData(polyData)
-# writer.Write()
+polyData = vtk.vtkPolyData()
+polyData.DeepCopy(actor2.GetMapper().GetInput())
+transform = vtk.vtkTransform()
+transform.SetMatrix(actor2.GetMatrix())
+fil = vtk.vtkTransformPolyDataFilter()
+fil.SetTransform(transform)
+fil.SetInputDataObject(polyData)
+fil.Update()
+polyData.DeepCopy(fil.GetOutput())
+
+writer = vtk.vtkPLYWriter()
+writer.SetFileTypeToASCII()
+writer.SetColorModeToDefault()
+filename = r'C:\Users\mariskawesseli\Documents\GitLab\femur_lig_ply_col2.ply'
+writer.SetFileName(filename)
+writer.SetInputData(polyData)
+writer.Write()
 
 # import pandas as pd
 # pd.DataFrame(color).to_clipboard()
\ No newline at end of file
diff --git a/LigamentStudy/VisualizeMeanSSM.ipynb b/LigamentStudy/VisualizeMeanSSM.ipynb
index 69fd983..b5e7158 100644
--- a/LigamentStudy/VisualizeMeanSSM.ipynb
+++ b/LigamentStudy/VisualizeMeanSSM.ipynb
@@ -21,7 +21,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 1,
    "id": "0a85e4e1-b480-4f37-a5de-4d2098d01e66",
    "metadata": {
     "tags": []
@@ -45,14 +45,13 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 2,
    "id": "9ad2e1fd-d4a4-48d2-8966-671b64ce093d",
    "metadata": {},
    "outputs": [],
    "source": [
     "import os\n",
     "import vtk\n",
-    "from numpy import random\n",
     "import trimesh\n",
     "import numpy as np\n",
     "import seaborn as sns\n",
@@ -69,7 +68,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 3,
    "id": "4efddfba-c04e-4023-8f6c-fdca0a5f25de",
    "metadata": {},
    "outputs": [],
@@ -120,7 +119,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 4,
    "id": "8f40bba1-dcd8-4d37-9682-d7d19e21bea0",
    "metadata": {},
    "outputs": [],
@@ -159,7 +158,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 5,
    "id": "77dfff23-beb1-413e-b73b-c20a9e00245f",
    "metadata": {},
    "outputs": [],
@@ -179,12 +178,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 6,
    "id": "e418af21-2c51-441b-8e2c-94d800a1fd73",
    "metadata": {},
    "outputs": [],
    "source": [
-    "path = os.path.join(r'C:\\Users\\mariskawesseli\\Documents\\GitLab\\knee_ssm\\OAI\\Output/' + segment + r'_bone\\new_bone\\shape_models')"
+    "path = os.path.join(r'./data/' + segment)"
    ]
   },
   {
@@ -197,7 +196,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 7,
    "id": "154780b1-e011-4a8c-8935-03900ef064d2",
    "metadata": {},
    "outputs": [],
@@ -312,14 +311,14 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 12,
+   "execution_count": 11,
    "id": "c98f226f-a28a-4e44-a342-5771f4984d63",
    "metadata": {},
    "outputs": [
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "21bc484f709e4c0091acadc04799fac6",
+       "model_id": "dd208c6057af4d6696067520b403e9ee",
        "version_major": 2,
        "version_minor": 0
       },
@@ -347,6 +346,63 @@
     "plotter.show()"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "id": "cd376163",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "aee45e50ed6746e69e7e02c8392f989a",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "ViewInteractiveWidget(height=768, layout=Layout(height='auto', width='100%'), width=1024)"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "spheres=[]\n",
+    "plotter = pv.Plotter()\n",
+    "for i in range(0,len(points_lig)):\n",
+    "    spheres.append(pv.Sphere(center=points_lig[i], radius=0.25))\n",
+    "    cols = np.tile(rgb_col[i], (spheres[i].number_of_points,1))\n",
+    "    spheres[i][\"colors\"] = cols\n",
+    "    plotter.add_mesh(spheres[i])\n",
+    "\n",
+    "plotter.add_actor(bone_actor)\n",
+    "plotter.add_actor(legend)\n",
+    "pv.set_plot_theme(\"document\")\n",
+    "plotter.show()  # show the two spheres from two PolyData\n",
+    "\n",
+    "plotter.export_html(r\"C:\\Users\\mariskawesseli\\Documents\\GitLab\\2022_JCWMSK_tutorials\\SSMfemur.html\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "id": "990d654f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# plotter = pv.Plotter(window_size=(900, 900),notebook=True)\n",
+    "# mesh= pv.read(r\"C:\\Users\\mariskawesseli\\Documents\\GitLab\\femur_lig_ply_col.ply\")\n",
+    "# scalars = mesh['RGBA']\n",
+    "# plotter.add_actor(bone_actor)\n",
+    "# plotter.add_mesh(mesh, show_scalar_bar=False, scalars=scalars[:,0:3])\n",
+    "# plotter.add_actor(legend)\n",
+    "# pv.set_plot_theme(\"document\")\n",
+    "# plotter.show()\n",
+    "\n",
+    "# plotter.export_html(r\"C:\\Users\\mariskawesseli\\Documents\\GitLab\\2022_JCWMSK_tutorials\\SSMfemur.html\")"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "9ae0ddf1-b716-4e73-ab28-06015b2a8bc7",
@@ -367,7 +423,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 13,
+   "execution_count": 15,
    "id": "dedd30e3-4598-4b88-9d0d-e8b4b6fe5f26",
    "metadata": {},
    "outputs": [],
@@ -387,12 +443,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 14,
+   "execution_count": 16,
    "id": "eeb6abea-1e7c-42a3-be1c-6a0d1730ff58",
    "metadata": {},
    "outputs": [],
    "source": [
-    "path = os.path.join(r'C:\\Users\\mariskawesseli\\Documents\\GitLab\\knee_ssm\\OAI\\Output/' + segment + r'_bone\\new_bone\\shape_models')"
+    "path = os.path.join(r'./data/' + segment)"
    ]
   },
   {
@@ -405,7 +461,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 15,
+   "execution_count": 17,
    "id": "2445f518-0319-4075-800b-d35bdabe434a",
    "metadata": {},
    "outputs": [],
@@ -431,7 +487,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 16,
+   "execution_count": 18,
    "id": "606f3104-8996-4529-b101-796ac53e00c2",
    "metadata": {},
    "outputs": [],
@@ -455,14 +511,14 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 20,
+   "execution_count": 19,
    "id": "cd8b342b-39bd-4bcb-b0f6-f25d6d3edc2e",
    "metadata": {},
    "outputs": [
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "f21bda3655db48ddae6deabeb32d58a3",
+       "model_id": "278f323878ef45eba11f5d340b212af3",
        "version_major": 2,
        "version_minor": 0
       },
@@ -484,13 +540,71 @@
     "plotter.add_actor(legend)\n",
     "\n",
     "pv.set_plot_theme(\"document\")\n",
+    "\n",
     "plotter.show()"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": 20,
+   "id": "944d007a",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "85366469fecb4683975db877749f97ed",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "ViewInteractiveWidget(height=768, layout=Layout(height='auto', width='100%'), width=1024)"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "spheres=[]\n",
+    "plotter = pv.Plotter()\n",
+    "for i in range(0,len(points_lig)):\n",
+    "    spheres.append(pv.Sphere(center=points_lig[i], radius=0.25))\n",
+    "    cols = np.tile(rgb_col[i], (spheres[i].number_of_points,1))\n",
+    "    spheres[i][\"colors\"] = cols\n",
+    "    plotter.add_mesh(spheres[i])\n",
+    "\n",
+    "plotter.add_actor(bone_actor)\n",
+    "plotter.add_actor(legend)\n",
+    "pv.set_plot_theme(\"document\")\n",
+    "plotter.show()  # show the two spheres from two PolyData\n",
+    "\n",
+    "plotter.export_html(r\"C:\\Users\\mariskawesseli\\Documents\\GitLab\\2022_JCWMSK_tutorials\\SSMtibia.html\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 21,
+   "id": "0b196960-88d0-471c-9bec-a7231a969c85",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# plotter = pv.Plotter(window_size=(900, 900),notebook=True)\n",
+    "# mesh= pv.read(r\"C:\\Users\\mariskawesseli\\Documents\\GitLab\\tibia_lig_ply_col.ply\")\n",
+    "# scalars = mesh['RGBA']\n",
+    "# plotter.add_actor(bone_actor)\n",
+    "# plotter.add_mesh(mesh, show_scalar_bar=False, scalars=scalars[:,0:3])\n",
+    "# plotter.add_actor(legend)\n",
+    "# pv.set_plot_theme(\"document\")\n",
+    "# plotter.show()\n",
+    "\n",
+    "# plotter.export_html(r\"C:\\Users\\mariskawesseli\\Documents\\GitLab\\2022_JCWMSK_tutorials\\SSMtibia.html\")"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "7a65e4e1-74df-4b42-a3b2-e865922dd018",
+   "id": "3c364012",
    "metadata": {},
    "outputs": [],
    "source": []
@@ -513,6 +627,37 @@
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
    "version": "3.9.12"
+  },
+  "latex_envs": {
+   "LaTeX_envs_menu_present": true,
+   "autoclose": false,
+   "autocomplete": true,
+   "bibliofile": "biblio.bib",
+   "cite_by": "apalike",
+   "current_citInitial": 1,
+   "eqLabelWithNumbers": true,
+   "eqNumInitial": 1,
+   "hotkeys": {
+    "equation": "Ctrl-E",
+    "itemize": "Ctrl-I"
+   },
+   "labels_anchors": false,
+   "latex_user_defs": false,
+   "report_style_numbering": false,
+   "user_envs_cfg": false
+  },
+  "toc": {
+   "base_numbering": 1,
+   "nav_menu": {},
+   "number_sections": true,
+   "sideBar": true,
+   "skip_h1_title": false,
+   "title_cell": "Table of Contents",
+   "title_sidebar": "Contents",
+   "toc_cell": false,
+   "toc_position": {},
+   "toc_section_display": true,
+   "toc_window_display": false
   }
  },
  "nbformat": 4,
diff --git a/LigamentStudy/testVisualizeSSM.py b/LigamentStudy/testVisualizeSSM.py
new file mode 100644
index 0000000..d3273f9
--- /dev/null
+++ b/LigamentStudy/testVisualizeSSM.py
@@ -0,0 +1,441 @@
+# import pyvista as pv
+# mesh= pv.read(r"C:\Users\mariskawesseli\Documents\GitLab\femur_lig_ply_col.ply")
+# mesh.plot()
+
+# import pyvista as pv
+# import numpy as np
+# # Re cast PolyData because file was not properly saved
+# bad = pv.read(r"C:\Users\mariskawesseli\Documents\GitLab\femur_lig_ply_col.ply")
+# bad.plot()
+# mesh = pv.PolyData(bad.points)
+# # Plot it
+# scalars = bad['RGBA']
+# # mesh.plot(scalars=scalars)
+# mesh.plot(scalars=scalars[:,0:3])
+# mesh.plot(scalars=scalars)
+# mesh.plot(scalars=scalars, rgba=True)
+
+import sys
+import os
+import vtk
+from numpy import random
+import trimesh
+import numpy as np
+import seaborn as sns
+
+
+class VtkPointCloud:
+    def __init__(self, zMin=-10.0, zMax=10.0, maxNumPoints=1e6):
+        self.maxNumPoints = maxNumPoints
+        self.vtkPolyData = vtk.vtkPolyData()
+        self.clearPoints()
+        mapper = vtk.vtkPolyDataMapper()
+        mapper.SetInputData(self.vtkPolyData)
+        mapper.SetColorModeToDefault()
+        mapper.SetScalarRange(zMin, zMax)
+        mapper.SetScalarVisibility(1)
+        self.vtkActor = vtk.vtkActor()
+        self.vtkActor.SetMapper(mapper)
+
+    def addPoint(self, point):
+        if (self.vtkPoints.GetNumberOfPoints() < self.maxNumPoints):
+            pointId = self.vtkPoints.InsertNextPoint(point[:])
+            self.vtkDepth.InsertNextValue(point[2])
+            self.vtkCells.InsertNextCell(1)
+            self.vtkCells.InsertCellPoint(pointId)
+        else:
+            r = random.randint(0, self.maxNumPoints)
+            self.vtkPoints.SetPoint(r, point[:])
+        self.vtkCells.Modified()
+        self.vtkPoints.Modified()
+        self.vtkDepth.Modified()
+
+    def clearPoints(self):
+        self.vtkPoints = vtk.vtkPoints()
+        self.vtkCells = vtk.vtkCellArray()
+        self.vtkDepth = vtk.vtkDoubleArray()
+        self.vtkDepth.SetName('DepthArray')
+        self.vtkPolyData.SetPoints(self.vtkPoints)
+        self.vtkPolyData.SetVerts(self.vtkCells)
+        self.vtkPolyData.GetPointData().SetScalars(self.vtkDepth)
+        self.vtkPolyData.GetPointData().SetActiveScalars('DepthArray')
+
+
+def load_data(data, pointCloud):
+    # data = genfromtxt(filename, dtype=float, usecols=[0, 1, 2])
+    for k in range(size(data, 0)):
+        point = data[k]  # 20*(random.rand(3)-0.5)
+        pointCloud.addPoint(point)
+
+    return pointCloud
+
+
+def load_stl(filename):
+    reader = vtk.vtkSTLReader()
+    reader.SetFileName(filename)
+
+    mapper = vtk.vtkPolyDataMapper()
+    if vtk.VTK_MAJOR_VERSION <= 5:
+        mapper.SetInput(reader.GetOutput())
+    else:
+        mapper.SetInputConnection(reader.GetOutputPort())
+
+    actor = vtk.vtkActor()
+    actor.SetMapper(mapper)
+
+    return actor
+
+
+def create_pointcloud_polydata(points, colors=None, seg=None):
+    """https://github.com/lmb-freiburg/demon
+    Creates a vtkPolyData object with the point cloud from numpy arrays
+
+    points: numpy.ndarray
+        pointcloud with shape (n,3)
+
+    colors: numpy.ndarray
+        uint8 array with colors for each point. shape is (n,3)
+
+    Returns vtkPolyData object
+    """
+    vpoints = vtk.vtkPoints()
+    vpoints.SetNumberOfPoints(points.shape[0])
+    for i in range(points.shape[0]):
+        vpoints.SetPoint(i, points[i])
+    # vpoints.SetMarkerStyle(vtk.vtkPlotPoints.CIRCLE)
+
+    vpoly = vtk.vtkPolyData()
+
+    appendFilter = vtk.vtkAppendPolyData()
+    for i in range(points.shape[0]):
+        sphereSource = vtk.vtkSphereSource()
+        # spheres.SetThetaResolution(1)
+        # spheres.SetPhiResolution(1)
+        sphereSource.SetRadius(1)
+        sphereSource.SetCenter(vpoints.GetPoint(i))
+        sphereSource.Update()
+
+        appendFilter.AddInputData(sphereSource.GetOutput())
+
+    # vpoly.SetPoints(vpoints)
+    rgb_col = []
+    if not colors is None:
+        if seg == 'femur':
+            max_val = 8
+            color[112:len(color)] = (color[112:len(color)] / max_val) * 10
+        vcolors = vtk.vtkUnsignedCharArray()
+        vcolors.SetNumberOfComponents(3)
+        vcolors.SetName("Colors")
+        vcolors.SetNumberOfTuples(points.shape[0])
+        rgb_col = []
+        for i in range(points.shape[0]):
+            c = sns.color_palette("viridis_r", n_colors=101, as_cmap=False)
+            vcolors.SetTuple3(i, c[int(colors[i] * 10)][0] * 255, c[int(colors[i] * 10)][1] * 255,
+                              c[int(colors[i] * 10)][2] * 255)
+            rgb_col.append(
+                [c[int(colors[i] * 10)][0] * 255, c[int(colors[i] * 10)][1] * 255, c[int(colors[i] * 10)][2] * 255])
+            # print(i, c[int(colors[i] - 1)][0], c[int(colors[i] - 1)][1], c[int(colors[i] - 1)][2])
+            # c = rgb(1,10,colors[i])
+            # vcolors.SetTuple3(i, c[0], c[1], c[2])
+        vpoly.GetPointData().SetScalars(vcolors)
+
+        actor.GetProperty().SetColor(color)
+
+    vcells = vtk.vtkCellArray()
+
+    for i in range(points.shape[0]):
+        vcells.InsertNextCell(1)
+        vcells.InsertCellPoint(i)
+
+    vpoly.SetVerts(vcells)
+
+    return vpoly, rgb_col
+
+
+def rgb(minimum, maximum, value):
+    minimum, maximum = float(minimum), float(maximum)
+    ratio = (value - minimum) / (maximum - minimum)  # 2 *
+    g = int(max(0, 255 * (1 - ratio)))
+    r = int(max(0, 255 * (ratio - 0)))
+    b = 0  # 255 - b - r
+    return r, g, b
+
+
+def createSpline(points):
+    vpoints = vtk.vtkPoints()
+    vpoints.SetNumberOfPoints(points.shape[0])
+    for i in range(points.shape[0]):
+        vpoints.SetPoint(i, points[i])
+
+    spline = vtk.vtkParametricSpline()
+    spline.SetPoints(vpoints)
+
+    functionSource = vtk.vtkParametricFunctionSource()
+    functionSource.SetParametricFunction(spline)
+    functionSource.Update()
+
+    # Create a mapper
+    mapper = vtk.vtkPolyDataMapper()
+    mapper.SetInputConnection(functionSource.GetOutputPort())
+
+    # Create an actor
+    actor = vtk.vtkActor()
+    actor.SetMapper(mapper)
+
+    return actor
+
+
+if __name__ == '__main__':
+    center_tibia = np.concatenate((np.arange(131), np.arange(470 - 341) + 341))  # PCL + ACL
+    center_femur = np.concatenate((np.arange(112), np.arange(341 - 263) + 263))  # PCL + ACL
+    # center_femur = np.concatenate((np.arange(64), np.arange(101 - 68) + 68))  # PCL + ACL
+    center_only = 1
+    subjects = [100]  # [100]  # ['9','13','19','23','26','29','32','35','37','41'] #, S0 [100]
+
+    segments = ['femur']  # 'femur',
+    ligaments_fem = [[1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
+                     [6, 5, 6, 6, 6, 6, 4, 4, 5, 5],
+                     [3, 2, 5, 3, 3, 2, 2, 0, 3, 3],
+                     [0, 8, 0, 0, 0, 0, 0, 0, 0, 0],  # MCLd2
+                     [7, 3, 7, 7, 7, 5, 7, 6, 7, 0],
+                     [0, 0, 8, 0, 0, 0, 0, 0, 0, 0],  # POL2
+                     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],  # POL3
+                     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],  # POL4
+                     [4, 6, 3, 5, 4, 0, 0, 3, 4, 4],
+                     [5, 7, 4, 4, 5, 7, 6, 5, 6, 6],
+                     [2, 4, 2, 2, 2, 3, 3, 2, 2, 2]]
+
+    ligaments_tib = [[5, 7, 6, 5, 3, 4, 4, 5, 5, 4],
+                     [1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
+                     [3, 3, 8, 3, 5, 3, 5, 0, 3, 3],
+                     [0, 4, 0, 0, 0, 0, 0, 0, 0, 0],  # MCLd2
+                     [4, 5, 3, 4, 4, 5, 3, 2, 4, 0],
+                     [0, 6, 4, 0, 0, 0, 0, 0, 0, 0],  # POL2
+                     [0, 0, 5, 0, 0, 0, 0, 0, 0, 0],  # POL3
+                     [0, 0, 7, 0, 0, 0, 0, 0, 0, 0],  # POL4
+                     [6, 8, 9, 6, 6, 6, 6, 6, 6, 5],
+                     [2, 2, 2, 2, 2, 2, 2, 3, 2, 2],
+                     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]
+
+    ligaments_fib = [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],  # PCL
+                     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],  # MCLp
+                     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],  # MCLd
+                     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],  # MCLd2
+                     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],  # POL
+                     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],  # POL2
+                     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],  # POL3
+                     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],  # POL4
+                     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],  # ACL
+                     [2, 2, 2, 2, 2, 2, 2, 3, 2, 2],  # LCL
+                     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]  # POP
+
+    for segment in segments:
+        SSMpoints = [[] for i in range(11)]
+        if segment == 'tibia':
+            center = center_tibia
+        elif segment == 'femur':
+            center = center_femur
+
+        for ind in range(0, 11):
+            SSMpoints[ind] = [[] for i in range(10)]
+
+        for ind, subject in enumerate(subjects):
+            if subject == 100:
+                path = os.path.join(
+                    r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone\shape_models')
+            elif subject == 'S0':
+                path = os.path.join(r'C:\Users\mariskawesseli\Documents\LigamentStudy\MRI\S0_prelim')
+            else:
+                path = os.path.join(r"C:\Users\mariskawesseli\Documents\LigamentStudy\ImageData", str(subject))
+
+            if subject in [9, 13, 26, 29, 32]:
+                side = 'R'
+                reflect = ''
+            else:
+                side = 'L'
+                reflect = '.reflect'
+
+            # points = trimesh.load_mesh(r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output\femur_bone\new_bone\shape_models\meanshape_bone_no_lig.xyz')
+            # point_cloud = create_pointcloud_polydata(points)
+            # pointCloud = VtkPointCloud()
+            # pointCloud = load_data(point_cloud, pointCloud)
+            # points_lig = trimesh.load_mesh(r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output\femur_bone\new_bone\shape_models\meanshape_ligs.xyz')
+            if subject == 100:
+                # points_lig = trimesh.load_mesh(path + '\meanshape_ligs.xyz')
+                # point_cloud_lig = create_pointcloud_polydata(points_lig)
+                points_lig = trimesh.load_mesh(path + '\meanshape_ligs_color.xyz')
+                color = np.loadtxt(path + r'\meanshape_ligs_color.xyz')[:, 3]
+
+                if center_only == 1:
+                    points_lig = points_lig[center]
+                    color = color[center]
+                point_cloud_lig, rgb_col = create_pointcloud_polydata(points_lig, colors=color, seg=segment)
+                bone_actor = load_stl(path + '/mean_shape.stl')
+                bone_actor.GetProperty().SetOpacity(1.0)
+
+                mesh = trimesh.load_mesh(path + '/mean_shape.stl')
+                # dist = trimesh.proximity.nearby_faces(mesh, np.squeeze(np.asarray(points_lig[np.argwhere(color >= 8)])))
+                dist3 = trimesh.proximity.closest_point_naive(mesh, np.squeeze(
+                    np.asarray(points_lig[np.argwhere(color >= 7)])), tol=1.0)
+
+                # faces = np.unique(np.asarray([item for sublist in dist for item in sublist]))
+                faces = np.unique(np.asarray([item for sublist in dist3[3] for item in sublist]))
+                mesh.update_faces(faces)
+                mesh.export(path + '/mean_shape_80percsurf.stl')
+                surf_actor = load_stl(path + '/mean_shape_80percsurf.stl')
+            else:
+                # points_lig = trimesh.load_mesh(path + '\SSM_' + segment + '_areas.xyz')  #_pred_points_color
+                # point_cloud_lig = create_pointcloud_polydata(points_lig)
+                points_lig = trimesh.load_mesh(
+                    path + '\SSM_' + segment + '_pred_points_color.xyz')  # _pred_points_color
+                color = np.loadtxt(path + '\SSM_' + segment + '_pred_points_color.xyz')[:,
+                        3]  # _areas _short_areas _pred_points
+                if center_only == 1:
+                    points_lig = points_lig[center]
+                    # color = color[center]
+                point_cloud_lig = create_pointcloud_polydata(points_lig, seg=segment)  # ,color colors=color,
+                if subject == 'S0':
+                    # bone_actor = load_stl(path + '/bone_femur2_2_bone_rot.stl')
+                    # bone_actor = load_stl(path + '/bone_tibia_2_bone_rot.stl')
+                    bone_actor = load_stl(path + '/bone_fibula_1_tissue_rot.stl')
+                else:
+                    bone_actor = load_stl(
+                        path + '/Segmentation_' + segment + '_resample.stl')  # '/SSM_' + segment + '_reconstruct_transform_icp.stl'
+                    if segment == 'fibula':
+                        segment_temp = 'tibia'
+                    else:
+                        segment_temp = segment
+                    # if center_only == 1:
+                    # 	wire_actor = load_stl(path + '/Segmentation_' + segment_temp + '_wires1.stl')
+                    # 	wire_actor2 = load_stl(path + '/Segmentation_' + segment_temp + '_wires3.stl')
+                    # 	wire_actor2.GetProperty().SetColor(1, 1, 0)
+                    # else:
+                    wire_actor = load_stl(path + '/Segmentation_' + segment_temp + '_wires.stl')
+                    wire_actor.GetProperty().SetColor(1, 1, 0)
+                bone_actor.GetProperty().SetOpacity(0.75)
+
+                points_bone = trimesh.load_mesh(path + '\SSM_' + segment + '_transform_icp.xyz')
+                point_cloud_bone = create_pointcloud_polydata(points_bone)
+
+                # orders = np.load(r'C:\Users\mariskawesseli\Documents\LigamentStudy\ImageData\occurances_order.npy')
+
+                mapper = vtk.vtkPolyDataMapper()
+                mapper.SetInputData(point_cloud_bone)
+                actor = vtk.vtkActor()
+                actor.SetMapper(mapper)
+                actor.GetProperty().SetColor(0, 0, 0)
+                actor.GetProperty().SetPointSize(2)
+            # actor.GetProperty().SetOpacity(1.0)
+
+            # spline_actor = createSpline(np.squeeze(np.asarray(points_lig[np.argwhere(color >= 8)])))
+            bone_actor.GetProperty().SetColor(0.89, 0.85, 0.79)
+            # bone_actor.GetProperty().LightingOff()
+            mapper2 = vtk.vtkPolyDataMapper()
+            mapper2.SetInputData(point_cloud_lig)
+            actor2 = vtk.vtkActor()
+            actor2.SetMapper(mapper2)
+            actor2.GetProperty().RenderPointsAsSpheresOn()
+            actor2.GetProperty().SetColor(1, 0, 0)
+            actor2.GetProperty().SetPointSize(7.5)
+
+            c = sns.color_palette("viridis_r", n_colors=101, as_cmap=False)
+            lut = vtk.vtkLookupTable()
+            lut.SetNumberOfColors(11)
+            lut.SetTableRange(1, 11)
+            for j in range(0, 11):
+                lut.SetTableValue(int(j * 1), c[j * 10][0], c[j * 10][1], c[j * 10][2])
+                # print(int(j*1), c[j*10-1][0], c[j*10-1][1], c[j*10-1][2])
+
+            j = 10 - 1
+            surf_col = [c[j][0], c[j][1], c[j][2]]
+            surf_col = [169 / 255, 169 / 255, 169 / 255]
+            surf_actor.GetProperty().SetColor(surf_col)
+            surf_actor.GetProperty().SetOpacity(1.0)
+
+            legend = vtk.vtkScalarBarActor()
+            legend.SetOrientationToHorizontal()
+            labelFormat = vtk.vtkTextProperty()
+            labelFormat.SetFontSize(16)
+            titleFormat = vtk.vtkTextProperty()
+            titleFormat.SetFontSize(8)
+            legend.SetLabelTextProperty(labelFormat)
+            # legend.SetTitleTextProperty(titleFormat)
+
+            legend.SetNumberOfLabels(11)
+            lut.SetTableRange(0, 100)
+            legend.SetLookupTable(lut)
+            # pos = legend.GetPositionCoordinate()
+            # pos.SetCoordinateSystemToNormalizedViewport()
+
+            legend.SetTitle("% of specimens \n")
+            legend.SetLabelFormat("%1.0f")
+            legend.SetUnconstrainedFontSize(1)
+
+            text_prop_cb = legend.GetLabelTextProperty()
+            text_prop_cb.SetFontFamilyAsString('Arial')
+            text_prop_cb.SetFontFamilyToArial()
+            text_prop_cb.SetColor(0, 0, 0)
+            # text_prop_cb.SetFontSize(500)
+            text_prop_cb.ShadowOff()
+            legend.SetLabelTextProperty(text_prop_cb)
+            # legend.SetMaximumWidthInPixels(75)
+            # legend.SetMaximumHeightInPixels(300)
+            legend.SetMaximumWidthInPixels(300)
+            legend.SetMaximumHeightInPixels(75)
+            legend.SetTitleTextProperty(text_prop_cb)
+            # legend.SetPosition(0.85,0.5)
+            legend.SetPosition(0.5, 0.85)
+
+            # Renderer
+            renderer = vtk.vtkRenderer()
+            # renderer.AddActor(actor)
+            renderer.AddActor(actor2)
+            renderer.AddActor(bone_actor)
+            # renderer.AddActor(spline_actor)
+            # renderer.AddActor(surf_actor)
+            if not subject == 100 and not subject == 'S0':
+                renderer.AddActor(wire_actor)
+                # renderer.AddActor(wire_actor2)
+            renderer.AddActor(legend)
+            # renderer.SetBackground(.2, .3, .4)
+            renderer.SetBackground(1.0, 1.0, 1.0)
+            renderer.ResetCamera()
+            # light = vtk.vtkLight()
+            # light.SetIntensity(1)
+            # renderer.AddLight(light)
+
+            # Render Window
+            renderWindow = vtk.vtkRenderWindow()
+            renderWindow.AddRenderer(renderer)
+            renderWindow.SetSize(750, 750)
+
+            # Interactor
+            renderWindowInteractor = vtk.vtkRenderWindowInteractor()
+            renderWindowInteractor.SetRenderWindow(renderWindow)
+            renderWindowInteractor.GetInteractorStyle().SetCurrentStyleToTrackballCamera()
+
+            # Begin Interaction
+            renderWindow.Render()
+            renderWindow.SetWindowName("XYZ Data Viewer " + str(subject))
+            renderWindowInteractor.Start()
+
+polyData = vtk.vtkPolyData()
+polyData.DeepCopy(actor2.GetMapper().GetInput())
+transform = vtk.vtkTransform()
+transform.SetMatrix(actor2.GetMatrix())
+fil = vtk.vtkTransformPolyDataFilter()
+fil.SetTransform(transform)
+fil.SetInputDataObject(polyData)
+fil.Update()
+polyData.DeepCopy(fil.GetOutput())
+
+writer = vtk.vtkPLYWriter()
+writer.SetFileTypeToASCII()
+writer.SetColorModeToDefault()
+filename = r'C:\Users\mariskawesseli\Documents\GitLab\femur_lig_ply_col2.ply'
+writer.SetFileName(filename)
+writer.SetInputData(polyData)
+writer.Write()
+
+# import pandas as pd
+# pd.DataFrame(color).to_clipboard()
\ No newline at end of file
-- 
GitLab