From d822ddc9cbbbe5a6a944050a7f92811cbf9b8bf7 Mon Sep 17 00:00:00 2001
From: Mariska Wesseling <M.G.H.Wesseling@tudelft.nl>
Date: Fri, 18 Mar 2022 15:08:06 +0100
Subject: [PATCH] Added Blumensaat line to ligament study. Updates to meshing
code.
---
LigamentStudy/BlumensaatLine.py | 286 +++++++++++++++++++
LigamentStudy/VisualiseSSM.py | 66 ++++-
LigamentStudy/VisualizeProjectedCentroids.py | 2 +-
LigamentStudy/remesh.py | 65 +++++
4 files changed, 408 insertions(+), 11 deletions(-)
create mode 100644 LigamentStudy/BlumensaatLine.py
create mode 100644 LigamentStudy/remesh.py
diff --git a/LigamentStudy/BlumensaatLine.py b/LigamentStudy/BlumensaatLine.py
new file mode 100644
index 0000000..c63da3c
--- /dev/null
+++ b/LigamentStudy/BlumensaatLine.py
@@ -0,0 +1,286 @@
+# Find most anterior edge of the femoral notch roof - representation Blumensaat line for 3D shapes
+# https://journals.lww.com/jbjsjournal/Fulltext/2010/06000/The_Location_of_Femoral_and_Tibial_Tunnels_in.10.aspx?__hstc=215929672.82af9c9a98fa600b1bb630f9cde2cb5f.1528502400314.1528502400315.1528502400316.1&__hssc=215929672.1.1528502400317&__hsfp=1773666937&casa_token=BT765BcrC3sAAAAA:Vu9rn-q5ng4c8339KQuq2mGZDgrAgBStwvn4lvYEbvCgvKQZkbJL24hWbKFdnHTc8VBmAIXA3HVvuWg22-9Mvwv1sw
+# https://www.dropbox.com/sh/l7pd43t7c4hrjdl/AABkncBbleifnpLDKSDDc0dCa/D3%20-%20Dimitriou%202020%20-%20Anterior%20cruciate%20ligament%20bundle%20insertions%20vary.pdf?dl=0
+
+import trimesh
+import numpy as np
+import os
+import math
+import pandas as pd
+import pymeshlab
+
+
+def findIntersection(x1, y1, x2, y2, x3, y3, x4, y4):
+ px = ((x1 * y2 - y1 * x2) * (x3 - x4) - (x1 - x2) * (x3 * y4 - y3 * x4)) / (
+ (x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4))
+ py = ((x1 * y2 - y1 * x2) * (y3 - y4) - (y1 - y2) * (x3 * y4 - y3 * x4)) / (
+ (x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4))
+
+ ang = math.atan2(py - y3, px - x3) - math.atan2(y1 - y3, x1 - x3)
+
+ l = math.cos(ang)*np.linalg.norm(np.asarray((x3,y3))-np.asarray((x4,y4)))
+
+ return l, px, py
+
+
+ligaments_fem = [[1,1,1,1,1,1,1,1,1,1], # PCL
+ [6,5,6,6,6,6,4,4,5,5], # MCLp
+ [3,2,5,3,3,2,2,0,3,3], # MCLd
+ [0,8,0,0,0,0,0,0,0,0], # MCLd2
+ [7,3,7,7,7,5,7,6,7,0], # POL
+ [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], # ACL
+ [5,7,4,4,5,7,6,5,6,6], # LCL
+ [2,4,2,2,2,3,3,2,2,2]] # POP
+
+ligaments = ligaments_fem
+
+# find most ant point in yz plane
+subjects = [100] # [9,13,19,23,26,29,32,35,37,41] #
+lig = 'PCL'
+segment = 'femur'
+
+d = []
+h = []
+h_centriods = []
+d_centriods = []
+for ind, subject in enumerate(subjects):
+ if subject in [9, 13, 26, 29, 32]:
+ side = 'R'
+ reflect = ''
+ else:
+ side = 'L'
+ reflect = '.reflect'
+
+ if subject == 100:
+ path = os.path.join(
+ r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone\shape_models\mean_shape_rot.stl')
+ side = 'R'
+ else:
+ path = os.path.join(r"C:\Users\mariskawesseli\Documents\LigamentStudy\ImageData", str(subject), 'Segmentation_femur_transform.STL')
+
+ mesh = trimesh.load_mesh(path)
+ verts = mesh.vertices
+
+ # posterior_mesh = trimesh.intersections.slice_mesh_plane(mesh, (0,-1,0), (0,20,0), cached_dots=None, return_both=False)
+ # posterior_mesh.show()
+
+ # find blumensaat line
+ lines, to_3D, face_index = trimesh.intersections.mesh_multiplane(mesh, (0,0,0), (1,0,0), heights=np.linspace(-10, 10, 21))
+ dist_point = []
+ prox_point = []
+ for i in range(0,len(face_index)):
+ plane_verts = np.unique(mesh.faces[face_index[i]])
+ plane_points = mesh.vertices[plane_verts]
+
+ goon = 1
+ tel = 2
+ while goon == 1:
+ min_z = np.where(plane_points[:,2] == np.partition(plane_points[:,2], tel)[tel])
+ y_min = plane_points[min_z][0][1]
+ min_z2 = np.where(plane_points[:,2] == np.partition(plane_points[:,2], tel+1)[tel+1])
+ y_min2 = plane_points[min_z2][0][1]
+ if y_min-y_min2 > -15:
+ goon = 1
+ tel += 1
+ else:
+ goon = 0
+ dist_point.append(plane_points[min_z][0])
+ min_y = np.where(plane_points[:,1] == plane_points[:,1].min())
+ prox_point.append(plane_points[min_y][0])
+
+ most_ant_ind1 = np.asarray(dist_point)[:, 1].argmax()
+ most_ant_ind2 = np.asarray(prox_point)[:, 1].argmax()
+
+ p1 = []
+ p2 = []
+ if most_ant_ind1 == most_ant_ind2:
+ p1.append(dist_point[most_ant_ind1])
+ p2.append(prox_point[most_ant_ind1])
+ print('equal')
+ else:
+ p1.append(dist_point[most_ant_ind2])
+ p2.append(prox_point[most_ant_ind2])
+ print('not equal')
+
+ if side == 'R':
+ if lig == 'ACL':
+ lateral_mesh = trimesh.intersections.slice_mesh_plane(mesh, (1,0,0), (0, 0, 0), cached_dots=None, return_both=False)
+ else:
+ lateral_mesh = trimesh.intersections.slice_mesh_plane(mesh, (-1, 0, 0), (0, 0, 0), cached_dots=None,
+ return_both=False)
+ else:
+ if lig == 'ACL':
+ lateral_mesh = trimesh.intersections.slice_mesh_plane(mesh, (-1, 0, 0), (0, 0, 0), cached_dots=None,
+ return_both=False)
+ else:
+ lateral_mesh = trimesh.intersections.slice_mesh_plane(mesh, (1, 0, 0), (0, 0, 0), cached_dots=None,
+ return_both=False)
+
+ # find height
+ vec1 = (p1[0][0] - p2[0][0], p1[0][1] - p2[0][1], p1[0][2] - p2[0][2])
+ norm = np.sqrt(vec1[0] ** 2 + vec1[1] ** 2 + vec1[2] ** 2)
+ direction = [vec1[0] / norm, vec1[1] / norm, vec1[2] / norm]
+
+
+ # segments = np.asarray([p1[-1], p2[-1]])
+ # p = trimesh.load_path(segments)
+
+ # trimesh.path.segments.parameters_to_segments(p1[-1], -1*direction, ((0,0,0),(0,1,0)))
+ # trimesh.path.segments.segments_to_parameters(np.asarray(segments))
+
+ # posterior_mesh = trimesh.intersections.slice_mesh_plane(mesh, direction, (0,0,10), cached_dots=None, return_both=False)
+
+ lines, to_3D, face_index = trimesh.intersections.mesh_multiplane(lateral_mesh, (0,0,0), direction, heights=np.linspace(-10, 10, 21))
+
+ dist = []
+ p3 = []
+ p4 = []
+ p1_2d = p1[-1][1:3]
+ p2_2d = p2[-1][1:3]
+ for i in range(0,len(face_index)):
+ plane_verts = np.unique(lateral_mesh.faces[face_index[i]])
+ plane_points = lateral_mesh.vertices[plane_verts]
+
+ min_y = np.where(plane_points[:,1] == np.partition(plane_points[:,1], 0)[0])
+ max_y = np.where(plane_points[:,1] == np.partition(plane_points[:,1], -1)[-1])
+
+ p3.append(plane_points[min_y][0])
+ p4.append(plane_points[max_y][0])
+ dist.append(np.linalg.norm(np.cross(p2_2d-p1_2d, p1_2d-p3[i][1:3]))/np.linalg.norm(p2_2d-p1_2d))
+ # dist.append(np.linalg.norm(plane_points[min_y][0]-plane_points[max_y][0]))
+
+ # segments = np.asarray([p3[np.asarray(dist).argmax()], p4[np.asarray(dist).argmax()]])
+ # p_dist = trimesh.load_path(segments)
+ dist1 = dist
+ p3_2d = p3[np.asarray(dist1).argmax()][1:3]
+ h.append(np.linalg.norm(np.cross(p2_2d-p1_2d, p1_2d-p3_2d))/np.linalg.norm(p2_2d-p1_2d))
+
+ # find depth
+ # lateral_mesh.show()
+ lines, to_3D, face_index = trimesh.intersections.mesh_multiplane(lateral_mesh, (0, 0, 0), direction,
+ heights=np.linspace(-30, -5, 41))
+
+ dist = []
+ p6 = []
+ p4 = []
+ p1_2d = p1[-1][1:3]
+ p2_2d = p2[-1][1:3]
+ for i in range(0, len(face_index)):
+ plane_verts = np.unique(lateral_mesh.faces[face_index[i]])
+ plane_points = lateral_mesh.vertices[plane_verts]
+
+ min_y = np.where(plane_points[:, 1] == np.partition(plane_points[:, 1], 0)[0])
+ max_y = np.where(plane_points[:, 1] == np.partition(plane_points[:, 1], -1)[-1])
+
+ p6.append(plane_points[min_y][0])
+ p4.append(plane_points[max_y][0])
+ dist.append(np.linalg.norm(np.cross(p2_2d - p1_2d, p1_2d - p6[i][1:3])) / np.linalg.norm(p2_2d - p1_2d))
+ # dist.append(np.linalg.norm(plane_points[min_y][0]-plane_points[max_y][0]))
+
+ jump_ind = np.where(np.diff(np.asarray(p6), axis=0)[:,1] == np.min(np.diff(np.asarray(p6), axis=0)[:,1]))[0][0]
+
+ # segments = np.asarray([p6[jump_ind+1], p4[jump_ind+1]])
+ # p_dist = trimesh.load_path(segments)
+
+ p6_2d = p6[jump_ind+1][1:3]
+ # min_z = lateral_mesh.vertices[np.argmin(lateral_mesh.vertices[:,2])]
+ # p5 = np.asarray(min_z)
+ # p5_2d = p5[1:3]
+
+ direction = np.asarray(direction) * -1
+ direction_perp = np.array((direction[0], -direction[2], direction[1]))
+
+ lines, to_3D, face_index = trimesh.intersections.mesh_multiplane(lateral_mesh, p1[0], direction_perp,
+ heights=np.linspace(0, 1, 1))
+ plane_verts = np.unique(lateral_mesh.faces[face_index[0]])
+ plane_points = lateral_mesh.vertices[plane_verts]
+ min_z = np.where(plane_points[:, 2] == np.partition(plane_points[:, 2], 0)[0])
+ p5 = plane_points[min_z][0]
+ p5_2d = p5[1:3]
+
+ l, px, py = findIntersection(p1_2d[0], p1_2d[1], p2_2d[0], p2_2d[1], p6_2d[0], p6_2d[1], p5_2d[0], p5_2d[1])
+ d.append(l)
+
+ # visualization
+ # p1[0][0] = 0
+ # p2[0][0] = 0
+ # p3[np.asarray(dist1).argmax()][0] = 0
+ # p4[jump_ind + 1][0] = 0
+ # p5[0] = 0
+ # p6[jump_ind + 1][0] = 0
+
+ points = trimesh.points.PointCloud(np.asarray((p1[0],p2[0],p6[jump_ind+1],p4[jump_ind+1],p5, p3[np.asarray(dist1).argmax()])), colors=None, metadata=None)
+ segments = np.asarray([p1[-1], p2[-1]])
+ p = trimesh.load_path(segments)
+ segments = np.asarray([p6[jump_ind+1], p5])
+ p_dist = trimesh.load_path(segments)
+
+ mesh.visual.face_colors[:] = np.array([227, 218, 201, 150])
+ mesh.visual.vertex_colors[:] = np.array([227, 218, 201, 150])
+ line = trimesh.path.segments.parameters_to_segments([p1[-1],p6[jump_ind+1],p3[np.asarray(dist1).argmax()],p5], [direction,direction_perp,direction,direction_perp],
+ np.array(((d[-1]-7,-7),(-10,h[-1]-10),(d[-1]-25,-25),(0,h[-1]))).astype(float))
+ trimesh.Scene([mesh, points, trimesh.load_path(np.squeeze(line))]).show()
+ # mesh.vertices[:, 0] = 0
+ # trimesh.Scene([mesh, points, trimesh.load_path(np.squeeze(line))]).show()
+
+# posterior_mesh = trimesh.intersections.slice_mesh_plane(mesh, direction, (0,-30,0), cached_dots=None, return_both=False)
+# posterior_mesh.show()
+ if subject == 100:
+ points_lig = trimesh.load_mesh(r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone\shape_models\meanshape_ligs_rot.xyz')
+ if lig == 'ACL':
+ center = np.arange(341 - 263) + 263 # ACL
+ else:
+ center = np.arange(112) # PCL
+ points_lig = points_lig[center]
+ # origin, xaxis, yaxis, zaxis = [0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]
+ # Rz = trimesh.transformations.rotation_matrix(180/np.pi, zaxis)
+ # points_lig.apply_transform(Rz)
+ for point in points_lig:
+ center_2d = point[1:3]
+ h_centriods.append(np.linalg.norm(np.cross(p2_2d - p1_2d, p1_2d - center_2d)) / np.linalg.norm(p2_2d - p1_2d))
+ l, px, py = findIntersection(p1_2d[0], p1_2d[1], p2_2d[0], p2_2d[1], center_2d[0], center_2d[1], p5_2d[0],
+ p5_2d[1])
+ d_centriods.append(l)
+ p_lig = trimesh.points.PointCloud(points_lig)
+ trimesh.Scene([mesh, points, p_lig, trimesh.load_path(np.squeeze(line))]).show()
+ else:
+ if lig == 'ACL':
+ lig_no = ligaments[8][ind]
+ elif lig == 'PCL':
+ lig_no = ligaments[0][ind]
+ if not lig_no == 0:
+ segment = 'femur'
+ path = os.path.join(r"C:\Users\mariskawesseli\Documents\LigamentStudy\ImageData", str(subject))
+
+ rot_mat = np.linalg.inv(np.loadtxt(path + '\Segmentation_' + segment + '_resample._ACS.txt'))
+ ms4 = pymeshlab.MeshSet()
+ ms4.load_new_mesh(path + '\Segmentation_' + segment + '_area' + str(lig_no) + '.stl')
+
+ ms4.apply_filter('flatten_visible_layers', deletelayer=True)
+ ms4.apply_filter('matrix_set_copy_transformation', transformmatrix=rot_mat)
+ geometric_measures = ms4.apply_filter('compute_geometric_measures')
+
+ # print('Surface area femur ligament' + str(lig_no) + ': ' + str(surface) + ' mm2')
+ center = geometric_measures['shell_barycenter']
+ center_2d = center[1:3]
+ h_centriods.append(np.linalg.norm(np.cross(p2_2d-p1_2d, p1_2d-center_2d))/np.linalg.norm(p2_2d-p1_2d))
+ l, px, py = findIntersection(p1_2d[0], p1_2d[1], p2_2d[0], p2_2d[1], center_2d[0], center_2d[1], p5_2d[0], p5_2d[1])
+ d_centriods.append(l)
+ else:
+ h_centriods.append(0)
+ d_centriods.append(0)
+
+[1-abs(i / j) for i, j in zip(d_centriods, d)]
+[i / j for i, j in zip(h_centriods, h)]
+
+d_centriods/np.asarray(d)
+h_centriods/np.asarray(h)
+
+np.mean(abs(np.asarray(d_centriods))/np.asarray(d))
+np.mean(h_centriods/np.asarray(h))
+
+
+
diff --git a/LigamentStudy/VisualiseSSM.py b/LigamentStudy/VisualiseSSM.py
index a5175ad..608bd61 100644
--- a/LigamentStudy/VisualiseSSM.py
+++ b/LigamentStudy/VisualiseSSM.py
@@ -93,10 +93,11 @@ def create_pointcloud_polydata(points, colors=None):
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=10, as_cmap=False)
vcolors.SetTuple3(i, c[int(colors[i] - 1)][0]*255, c[int(colors[i] - 1)][1]*255, c[int(colors[i] - 1)][2]*255)
+ rgb_col.append([c[int(colors[i] - 1)][0] * 255, c[int(colors[i] - 1)][1] * 255, c[int(colors[i] - 1)][2] * 255])
# print(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])
@@ -111,7 +112,7 @@ def create_pointcloud_polydata(points, colors=None):
vpoly.SetVerts(vcells)
- return vpoly
+ return vpoly, rgb_col
def rgb(minimum, maximum, value):
@@ -124,9 +125,13 @@ def rgb(minimum, maximum, value):
if __name__ == '__main__':
- subjects = ['S0'] #['9','13','19','23','26','29','32','35','37','41'] #, S0 [100] #
+ 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] # ['9','13','19','23','26','29','32','35','37','41'] #, S0 [100]
- segments = ['fibula'] #'femur',
+ segments = ['tibia'] #'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],
@@ -165,6 +170,11 @@ if __name__ == '__main__':
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)]
@@ -193,25 +203,38 @@ if __name__ == '__main__':
# point_cloud_lig = create_pointcloud_polydata(points_lig)
points_lig = trimesh.load_mesh(path + '\meanshape_ligs_color.xyz')
color = np.loadtxt(path + '\meanshape_ligs_color.xyz')[:, 3]
- point_cloud_lig = create_pointcloud_polydata(points_lig, color)
+
+ if center_only == 1:
+ points_lig = points_lig[center]
+ color = color[center]
+ point_cloud_lig, rgb_col = create_pointcloud_polydata(points_lig, color)
bone_actor = load_stl(path + '/mean_shape.stl')
+ bone_actor.GetProperty().SetOpacity(0.75)
# bone_actor.GetProperty().SetOpacity(0.75)
else:
- # points_lig = trimesh.load_mesh(path + '\SSM_' + segment + '_short_areas.xyz') #_pred_points_color
+ # 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 + '_short_pred_points.xyz') # _pred_points_color
- # color = np.loadtxt(path + '\SSM_' + segment + '_pred_points_color.xyz')[:,3] #_areas _short_areas _pred_points
- point_cloud_lig = create_pointcloud_polydata(points_lig) #,color
+ 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,color) #,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')
+ 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)
@@ -228,6 +251,7 @@ if __name__ == '__main__':
actor.GetProperty().SetColor(0,0,0)
actor.GetProperty().SetPointSize(2)
# actor.GetProperty().SetOpacity(1.0)
+
bone_actor.GetProperty().SetColor(0.89, 0.85, 0.79)
# bone_actor.GetProperty().LightingOff()
mapper2 = vtk.vtkPolyDataMapper()
@@ -272,6 +296,7 @@ if __name__ == '__main__':
renderer.AddActor(bone_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)
@@ -295,3 +320,24 @@ if __name__ == '__main__':
renderWindow.SetWindowName("XYZ Data Viewer " + str(subject))
renderWindowInteractor.Start()
+
+# polyData = vtk.vtkPolyData()
+# polyData.DeepCopy(bone_actor.GetMapper().GetInput())
+# transform = vtk.vtkTransform()
+# transform.SetMatrix(bone_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\femur_bone_ply.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/VisualizeProjectedCentroids.py b/LigamentStudy/VisualizeProjectedCentroids.py
index 99d8876..28aa072 100644
--- a/LigamentStudy/VisualizeProjectedCentroids.py
+++ b/LigamentStudy/VisualizeProjectedCentroids.py
@@ -113,7 +113,7 @@ color = ((0.75,1,0.5),
(1,0.5,0))
if __name__ == '__main__':
- subjects = [9] #9,13,19,23,26,29,32,35,37,41
+ subjects = [35] #9,13,19,23,26,29,32,35,37,41
segments = ['femur'] #'femur',
ligaments_fem = [[1,1,1,1,1,1,1,1,1,1],
diff --git a/LigamentStudy/remesh.py b/LigamentStudy/remesh.py
new file mode 100644
index 0000000..4ba854c
--- /dev/null
+++ b/LigamentStudy/remesh.py
@@ -0,0 +1,65 @@
+import trimesh
+import pymeshlab
+import os
+import numpy as np
+
+# mesh = trimesh.load_mesh(r'C:\Users\mariskawesseli\Documents\Data\OAI\segmentation\2019_ATEZ_MEDIA-Supplementary-Material-OAI-ZIB\OAI-ZIB\segmentation\segmentation_meshes\femur_cartilage\mesh\9005075.segmentation_masks_femoral_cartilage_R.ply')
+# trimesh.remesh.subdivide_to_size(mesh.vertices, mesh.faces, 5, max_iter=10, return_index=False)
+# trimesh.exchange.export.export_mesh(mesh,r'C:\Users\mariskawesseli\Documents\Data\OAI\segmentation\2019_ATEZ_MEDIA-Supplementary-Material-OAI-ZIB\OAI-ZIB\segmentation\segmentation_meshes\femur_cartilage\mesh_resample\9005075.segmentation_masks_femoral_cartilage_R.ply', 'ply')
+
+inputDir = r'C:\Users\mariskawesseli\Documents\Data\OAI\segmentation\2019_ATEZ_MEDIA-Supplementary-Material-OAI-ZIB\OAI-ZIB\segmentation'
+datasetName = "tibia_cartilage"
+mesh_dir = inputDir + r'/segmentation_meshes/' + datasetName + '/mesh/med/'
+mesh_dir_out = inputDir + r'/segmentation_meshes/' + datasetName + '/mesh_resample/med/'
+
+files_mesh = []
+for file in sorted(os.listdir(mesh_dir)):
+ files_mesh.append(mesh_dir + file)
+
+pt_to_use = r'C:\Users\mariskawesseli\Documents\Data\OAI\segmentation\2019_ATEZ_MEDIA-Supplementary-Material-OAI-ZIB\healthyKL_pts.txt'
+with open(pt_to_use) as f:
+ pts = f.readlines()
+pts = [i.split('\n')[0] for i in pts]
+pts_use = pts
+
+matches_mesh = []
+for pt in pts_use:
+ if any(pt in s for s in files_mesh):
+ matches_mesh.append([match for match in files_mesh if pt in match])
+
+files_mesh = [item for sublist in matches_mesh for item in sublist]
+
+#
+# for file in files_mesh:
+# ms6 = pymeshlab.MeshSet()
+# ms6.load_new_mesh(file)
+# ms6.apply_filter('uniform_mesh_resampling', cellsize=1.05, offset=0.5, mergeclosevert=False)
+# ms6.save_current_mesh(file_name=mesh_dir_out + file.split('/')[-1], save_textures=False)
+
+ # for femur
+ # m = ms6.current_mesh()
+ # fm = m.face_matrix()
+ #
+ # # ms6 = pymeshlab.MeshSet()
+ # # ms6.load_new_mesh(mesh_dir_out + file.split('/')[-1])
+ # ms6.apply_filter('simplification_quadric_edge_collapse_decimation',targetfacenum=np.max(fm),targetperc=0,qualitythr=0.3,preserveboundary=True,preservenormal=True,preservetopology=True)
+ # ms6.save_current_mesh(file_name=mesh_dir_out + file.split('/')[-1], save_textures=False)
+
+datasetName = "tibia_cartilage"
+side = 'med' # 'lat' #
+outputDirectory = r'C:/Users/mariskawesseli/Documents/GitLab/knee_ssm/OAI/Output/tibia_cartilage_' + side + '/groomed/'
+
+origin, xaxis, yaxis, zaxis = [0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]
+Rx = trimesh.transformations.rotation_matrix(np.pi/4, xaxis)
+
+mesh_reg = trimesh.load_mesh(outputDirectory + 'reference.ply')
+for file in files_mesh:
+ mesh1 = trimesh.load_mesh(mesh_dir_out + file.split('/')[-1])
+
+ T, cost = trimesh.registration.mesh_other(mesh1, mesh_reg, samples=500, scale=False, icp_first=10, icp_final=50)
+ mesh1.apply_transform(T)
+ mesh1.apply_transform(Rx)
+
+ mesh1.export(outputDirectory + 'meshes/' + file.split('/')[-1])
+
+
--
GitLab