From 016da880ccdf26bcab593ebee966ddafe7d75199 Mon Sep 17 00:00:00 2001
From: Mariska Wesseling <M.G.H.Wesseling@tudelft.nl>
Date: Thu, 18 Nov 2021 11:51:52 +0100
Subject: [PATCH] updates ligamentstudy
---
LigamentStudy/Analyses.py | 231 +++++++++++--------
LigamentStudy/AreaTest.py | 102 ++++++++
LigamentStudy/AttechmentArea.py | 206 +++++++++--------
LigamentStudy/HausdorffDistance.py | 136 +++++++++++
LigamentStudy/OAIdownload.py | 119 ++++++++++
LigamentStudy/ProjectCentroids.py | 213 +++++++++++++++++
LigamentStudy/SlicerExportXray.py | 173 ++++++++++++++
LigamentStudy/SlicerPositionBeam.py | 36 +++
LigamentStudy/VisualiseSSM.py | 105 ++++++---
LigamentStudy/VisualizeProjectedCentroids.py | 207 +++++++++++++++++
LigamentStudy/Xray.py | 140 ++++++-----
LigamentStudy/fitSSM.py | 220 +++++++++++-------
LigamentStudy/fitSSM_mri.py | 219 ++++++++++++++++++
LigamentStudy/scaleOsim.py | 170 ++++++++++++++
14 files changed, 1905 insertions(+), 372 deletions(-)
create mode 100644 LigamentStudy/AreaTest.py
create mode 100644 LigamentStudy/HausdorffDistance.py
create mode 100644 LigamentStudy/OAIdownload.py
create mode 100644 LigamentStudy/ProjectCentroids.py
create mode 100644 LigamentStudy/SlicerExportXray.py
create mode 100644 LigamentStudy/SlicerPositionBeam.py
create mode 100644 LigamentStudy/VisualizeProjectedCentroids.py
create mode 100644 LigamentStudy/fitSSM_mri.py
create mode 100644 LigamentStudy/scaleOsim.py
diff --git a/LigamentStudy/Analyses.py b/LigamentStudy/Analyses.py
index f1d74c4..913afa8 100644
--- a/LigamentStudy/Analyses.py
+++ b/LigamentStudy/Analyses.py
@@ -11,51 +11,58 @@ from openpyxl import load_workbook
# femur
# PCL: [1,1,1,1,1,1,1,1,1,1]
-# MCL-p: [6,5,6,6,6,6,4,4,5,5]
-# MCL-d: [3,2,5,3,3,2,2,-,3,3]
-# posterior oblique: [7,8,7,7,7,5,7,6,7,-]
-# ACL: [4,6,3,5,4,(4),(5),3,4,4]
+# MCL-superficial: [6,5,6,6,6,6,4,4,5,5]
+# MCL-deep: [3,2+8,5,3,3,2,2,-,3,3]
+# posterior oblique: [7,3,7+8,7,7,5,7,6,7,-]
+# ACL: [4,6,3,5,4,-(4),-(5),3,4,4]
# LCL (prox): [5,7,4,4,5,7,6,5,6,6]
# popliteus (dist): [2,4,2,2,2,3,3,2,2,2]
-# other med: [-,3,8,-,-,-,-,-,-,-]
# tibia
# PCL: [5,7,6,5,3,4,4,5,5,4]
-# MCL-p: [3,3,7,3,5,3,5,4,3,3]
-# MCL-d: [1,1,1,1,1,1,1,1,1,1]
-# posterior oblique: [4,5,3,4,4,5,3,2,4,3]
+# MCL-superficial: [1,1,1,1,1,1,1,1,1,1]
+# MCL-deep: [3,3+4,8,3,5,3,5,-(4),3,3]
+# posterior oblique: [4,5+6,3+4+5+7,4,4,5,3,2,4,-]
# ACL: [6,8,9,6,6,6,6,6,6,5]
# LCL: [2,2,2,2,2,2,2,3,2,2]
-# popliteus (dist): [-,-,-,-,-,-,-,-,-,-]
+# popliteus: [-,-,-,-,-,-,-,-,-,-]
subjects = [9,13,19,23,26,29,32,35,37,41] #9,13,19,23,26,29,32,35,41
-segments = ['femur','tibia'] #
+segments = ['tibia'] #'femur',
+short = 1
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],
- [7,8,7,7,7,5,7,6,7,0],
- [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],
- [0,3,8,0,0,0,0,0,0,0]]
+ [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],
- [3,3,7,3,5,3,5,4,3,3],
- [1,1,1,1,1,1,1,1,1,1],
- [4,5,3,4,4,5,3,2,4,3],
- [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],
- [0,0,0,0,0,0,0,0,0,0]]
+ [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]]
book = load_workbook(os.path.join(r"C:\Users\mariskawesseli\Documents\LigamentStudy\ImageData","surfaces.xlsx"))
writer = pd.ExcelWriter(os.path.join(r"C:\Users\mariskawesseli\Documents\LigamentStudy\ImageData","surfaces.xlsx"), engine='openpyxl')
writer.book = book
for segment in segments:
- surface = np.empty((8, 10))
+ surface = np.empty((11, 10))
surface[:] = np.nan
- center = np.empty((8, 10))
+ center = np.empty((11, 10))
center[:] = np.nan
ML_size = np.empty((1, 10))
ML_size[:] = np.nan
@@ -87,15 +94,15 @@ for segment in segments:
bb_max_lat[:] = np.nan
bb_min_lat = np.empty((3, 10))
bb_min_lat[:] = np.nan
- dist_to_edge = np.empty((8, 10, 3))
+ dist_to_edge = np.empty((11, 10, 3))
dist_to_edge[:] = np.nan
- perc_of_len = np.empty((8, 10, 3))
+ perc_of_len = np.empty((11, 10, 3))
perc_of_len[:] = np.nan
- perc_of_len_med = np.empty((8, 10, 3))
+ perc_of_len_med = np.empty((11, 10, 3))
perc_of_len_med[:] = np.nan
- perc_of_len_lat = np.empty((8, 10, 3))
+ perc_of_len_lat = np.empty((11, 10, 3))
perc_of_len_lat[:] = np.nan
- center = np.empty((8, 10, 3))
+ center = np.empty((11, 10, 3))
center[:] = np.nan
if segment == 'femur':
ligaments = ligaments_fem
@@ -156,7 +163,7 @@ for segment in segments:
bb_min_lat[:, ind] = np.min(ms5.current_mesh().vertex_matrix(), 0)
# determine surface area attachments
- for lig in range(0, 8):
+ for lig in range(0, 11):
lig_no = ligaments[lig][ind]
if not lig_no == 0:
ms4 = pymeshlab.MeshSet()
@@ -175,29 +182,38 @@ for segment in segments:
AP_size_lat[0, ind], AP_size_lat[0, ind], AP_size_lat[0, ind])
df = pd.DataFrame({'PCLx': perc_of_len[0,:,0],
- 'MCL-px': perc_of_len[1, :, 0],
- 'MCL-dx:': perc_of_len[2, :, 0],
- 'posterior obliquex': perc_of_len[3, :, 0],
- 'ACLx': perc_of_len[4, :, 0],
- 'LCL (prox)x': perc_of_len[5, :, 0],
- 'popliteus (dist)x': perc_of_len[6, :, 0],
- 'other medx': perc_of_len[7, :, 0],
+ 'MCL-sx': perc_of_len[1, :, 0],
+ 'MCL-d1x:': perc_of_len[2, :, 0],
+ 'MCL-d2x:': perc_of_len[3, :, 0],
+ 'posterior oblique1x': perc_of_len[4, :, 0],
+ 'posterior oblique2x': perc_of_len[5, :, 0],
+ 'posterior oblique3x': perc_of_len[6, :, 0],
+ 'posterior oblique4x': perc_of_len[7, :, 0],
+ 'ACLx': perc_of_len[8, :, 0],
+ 'LCLx': perc_of_len[9, :, 0],
+ 'popliteusx': perc_of_len[10, :, 0],
'PCLy': perc_of_len[0,:,1],
- 'MCL-py': perc_of_len[1, :, 1],
- 'MCL-dy:': perc_of_len[2, :, 1],
- 'posterior obliquey': perc_of_len[3, :, 1],
- 'ACLy': perc_of_len[4, :, 1],
- 'LCL (prox)y': perc_of_len[5, :, 1],
- 'popliteus (dist)y': perc_of_len[6, :, 1],
- 'other medy': perc_of_len[7, :, 1],
+ 'MCL-sy': perc_of_len[1, :, 1],
+ 'MCL-d1y:': perc_of_len[2, :, 1],
+ 'MCL-d2y:': perc_of_len[3, :, 1],
+ 'posterior oblique1y': perc_of_len[4, :, 1],
+ 'posterior oblique2y': perc_of_len[5, :, 1],
+ 'posterior oblique3y': perc_of_len[6, :, 1],
+ 'posterior oblique4y': perc_of_len[7, :, 1],
+ 'ACLy': perc_of_len[8, :, 1],
+ 'LCLy': perc_of_len[9, :, 1],
+ 'popliteusy': perc_of_len[10, :, 1],
'PCLz': perc_of_len[0,:,2],
- 'MCL-pz': perc_of_len[1,:,2],
- 'MCL-dz:': perc_of_len[2,:,2],
- 'posterior obliquez': perc_of_len[3,:,2],
- 'ACLz': perc_of_len[4,:,2],
- 'LCL (prox)z': perc_of_len[5,:,2],
- 'popliteus (dist)z': perc_of_len[6,:,2],
- 'other medz': perc_of_len[7,:,2]
+ 'MCL-sz': perc_of_len[1,:,2],
+ 'MCL-d1z:': perc_of_len[2,:,2],
+ 'MCL-d2z:': perc_of_len[3, :, 2],
+ 'posterior oblique1z': perc_of_len[4,:,2],
+ 'posterior oblique2z': perc_of_len[5, :, 2],
+ 'posterior oblique3z': perc_of_len[6, :, 2],
+ 'posterior oblique4z': perc_of_len[7, :, 2],
+ 'ACLz': perc_of_len[8,:,2],
+ 'LCLz': perc_of_len[9,:,2],
+ 'popliteusz': perc_of_len[10,:,2]
})
means = df.mean(skipna=True)
std = df.std(skipna=True)
@@ -211,21 +227,27 @@ for segment in segments:
summary_ave_data.to_excel(writer, sheet_name='perc_of_len ' + segment)
df = pd.DataFrame({'PCLy': perc_of_len_med[0, :, 1],
- 'MCL-py': perc_of_len_med[1, :, 1],
- 'MCL-dy:': perc_of_len_med[2, :, 1],
- 'posterior obliquey': perc_of_len_med[3, :, 1],
- 'ACLy': perc_of_len_lat[4, :, 1],
- 'LCL (prox)y': perc_of_len_lat[5, :, 1],
- 'popliteus (dist)y': perc_of_len_lat[6, :, 1],
- 'other medy': perc_of_len_med[7, :, 1],
+ 'MCL-sy': perc_of_len_med[1, :, 1],
+ 'MCL-d1y:': perc_of_len_med[2, :, 1],
+ 'MCL-d2y:': perc_of_len_med[3, :, 1],
+ 'posterior oblique1y': perc_of_len_med[4, :, 1],
+ 'posterior oblique2y': perc_of_len_med[5, :, 1],
+ 'posterior oblique3y': perc_of_len_med[6, :, 1],
+ 'posterior oblique4y': perc_of_len_med[7, :, 1],
+ 'ACLy': perc_of_len_lat[8, :, 1],
+ 'LCLy': perc_of_len_lat[9, :, 1],
+ 'popliteusy': perc_of_len_lat[10, :, 1],
'PCLz': perc_of_len_med[0, :, 2],
- 'MCL-pz': perc_of_len_med[1, :, 2],
- 'MCL-dz:': perc_of_len_med[2, :, 2],
- 'posterior obliquez': perc_of_len_med[3, :, 2],
- 'ACLz': perc_of_len_lat[4, :, 2],
- 'LCL (prox)z': perc_of_len_lat[5, :, 2],
- 'popliteus (dist)z': perc_of_len_lat[6, :, 2],
- 'other medz': perc_of_len_med[7, :, 2]
+ 'MCL-sz': perc_of_len_med[1, :, 2],
+ 'MCL-d1z:': perc_of_len_med[2, :, 2],
+ 'MCL-d2z:': perc_of_len_med[3, :, 2],
+ 'posterior oblique1z': perc_of_len_med[4, :, 2],
+ 'posterior oblique2z': perc_of_len_med[5, :, 2],
+ 'posterior oblique3z': perc_of_len_med[6, :, 2],
+ 'posterior oblique4z': perc_of_len_med[7, :, 2],
+ 'ACLz': perc_of_len_lat[8, :, 2],
+ 'LCLz': perc_of_len_lat[9, :, 2],
+ 'popliteusz': perc_of_len_lat[10, :, 2]
})
means = df.mean(skipna=True)
std = df.std(skipna=True)
@@ -239,29 +261,38 @@ for segment in segments:
summary_ave_data.to_excel(writer, sheet_name='perc_of_len med_lat ' + segment)
df = pd.DataFrame({'PCLx': dist_to_edge[0,:,0],
- 'MCL-px': dist_to_edge[1, :, 0],
- 'MCL-dx:': dist_to_edge[2, :, 0],
- 'posterior obliquex': dist_to_edge[3, :, 0],
- 'ACLx': dist_to_edge[4, :, 0],
- 'LCL (prox)x': dist_to_edge[5, :, 0],
- 'popliteus (dist)x': dist_to_edge[6, :, 0],
- 'other medx': dist_to_edge[7, :, 0],
+ 'MCL-sx': dist_to_edge[1, :, 0],
+ 'MCL-d1x:': dist_to_edge[2, :, 0],
+ 'MCL-d2x:': dist_to_edge[3, :, 0],
+ 'posterior oblique1x': dist_to_edge[4, :, 0],
+ 'posterior oblique2x': dist_to_edge[5, :, 0],
+ 'posterior oblique3x': dist_to_edge[6, :, 0],
+ 'posterior oblique4x': dist_to_edge[7, :, 0],
+ 'ACLx': dist_to_edge[8, :, 0],
+ 'LCLx': dist_to_edge[9, :, 0],
+ 'popliteusx': dist_to_edge[10, :, 0],
'PCLy': dist_to_edge[0,:,1],
- 'MCL-py': dist_to_edge[1, :, 1],
- 'MCL-dy:': dist_to_edge[2, :, 1],
- 'posterior obliquey': dist_to_edge[3, :, 1],
- 'ACLy': dist_to_edge[4, :, 1],
- 'LCL (prox)y': dist_to_edge[5, :, 1],
- 'popliteus (dist)y': dist_to_edge[6, :, 1],
- 'other medy': dist_to_edge[7, :, 1],
+ 'MCL-sy': dist_to_edge[1, :, 1],
+ 'MCL-d1y:': dist_to_edge[2, :, 1],
+ 'MCL-d2y:': dist_to_edge[3, :, 1],
+ 'posterior oblique1y': dist_to_edge[4, :, 1],
+ 'posterior oblique2y': dist_to_edge[5, :, 1],
+ 'posterior oblique3y': dist_to_edge[6, :, 1],
+ 'posterior oblique4y': dist_to_edge[7, :, 1],
+ 'ACLy': dist_to_edge[8, :, 1],
+ 'LCLy': dist_to_edge[9, :, 1],
+ 'popliteusy': dist_to_edge[10, :, 1],
'PCLz': dist_to_edge[0,:,2],
- 'MCL-pz': dist_to_edge[1,:,2],
- 'MCL-dz:': dist_to_edge[2,:,2],
- 'posterior obliquez': dist_to_edge[3,:,2],
- 'ACLz': dist_to_edge[4,:,2],
- 'LCL (prox)z': dist_to_edge[5,:,2],
- 'popliteus (dist)z': dist_to_edge[6,:,2],
- 'other medz': dist_to_edge[7,:,2]
+ 'MCL-sz': dist_to_edge[1,:,2],
+ 'MCL-d1z:': dist_to_edge[2,:,2],
+ 'MCL-d2z:': dist_to_edge[3, :, 2],
+ 'posterior oblique1z': dist_to_edge[4,:,2],
+ 'posterior oblique2z': dist_to_edge[5, :, 2],
+ 'posterior oblique3z': dist_to_edge[6, :, 2],
+ 'posterior oblique4z': dist_to_edge[7, :, 2],
+ 'ACLz': dist_to_edge[8,:,2],
+ 'LCLz': dist_to_edge[9,:,2],
+ 'popliteusz': dist_to_edge[10,:,2]
})
means = df.mean(skipna=True)
std = df.std(skipna=True)
@@ -273,15 +304,23 @@ for segment in segments:
summary_ave_data = summary_ave_data.rename({10: 'mean', 11: 'std'}, axis='index')
summary_ave_data.to_excel(writer, sheet_name='distance_to_edge ' + segment)
-
+ MCLd = np.nansum([surface[2, :], surface[3, :]], 0)
+ MCLd[MCLd == 0] = 'nan'
+ pol = np.nansum([surface[4, :],surface[5, :],surface[6, :],surface[7, :]],0)
+ pol[pol == 0] = 'nan'
df = pd.DataFrame({'PCL': surface[0,:],
- 'MCL-p': surface[1,:],
- 'MCL-d:': surface[2,:],
- 'posterior oblique': surface[3,:],
- 'ACL': surface[4,:],
- 'LCL (prox)': surface[5,:],
- 'popliteus (dist)': surface[6,:],
- 'other med': surface[7,:]
+ 'MCL-s': surface[1,:],
+ 'MCL-d1:': surface[2,:],
+ 'MCL-d2:': surface[3, :],
+ 'MCL-d:': MCLd,
+ 'posterior oblique1': surface[4,:],
+ 'posterior oblique2': surface[5, :],
+ 'posterior oblique3': surface[6, :],
+ 'posterior oblique4': surface[7, :],
+ 'posterior oblique': pol,
+ 'ACL': surface[8,:],
+ 'LCL': surface[9,:],
+ 'popliteus': surface[10,:],
})
means = df.mean(skipna=True)
std = df.std(skipna=True)
diff --git a/LigamentStudy/AreaTest.py b/LigamentStudy/AreaTest.py
new file mode 100644
index 0000000..d940e83
--- /dev/null
+++ b/LigamentStudy/AreaTest.py
@@ -0,0 +1,102 @@
+import pandas as pd
+import os
+import trimesh
+import numpy as np
+import matplotlib.path as plt
+
+
+def heron(a, b, c):
+ s = (a + b + c) / 2
+ area = (s * (s - a) * (s - b) * (s - c)) ** 0.5
+ return area
+
+
+def distance3d(x1, y1, z1, x2, y2, z2):
+ a = (x1 - x2) ** 2 + (y1 - y2) ** 2 + (z1 - z2) ** 2
+ d = a ** 0.5
+ return d
+
+
+def area(x1, y1, z1, x2, y2, z2, x3, y3, z3):
+ a = distance3d(x1, y1, z1, x2, y2, z2)
+ b = distance3d(x2, y2, z2, x3, y3, z3)
+ c = distance3d(x3, y3, z3, x1, y1, z1)
+ A = heron(a, b, c)
+ return A
+
+
+# print("area of triangle is %r " %A)
+
+# A utility function to calculate area
+# of triangle formed by (x1, y1),
+# (x2, y2) and (x3, y3)
+
+# def area(x1, y1, x2, y2, x3, y3):
+# return abs((x1 * (y2 - y3) + x2 * (y3 - y1)
+# + x3 * (y1 - y2)) / 2.0)
+
+
+# A function to check whether point P(x, y)
+# lies inside the triangle formed by
+# A(x1, y1), B(x2, y2) and C(x3, y3)
+def isInside(p1, p2, p3, p):
+ x1 = p1[0]
+ y1 = p1[1]
+ z1 = p1[2]
+ x2 = p2[0]
+ y2 = p2[1]
+ z2 = p2[2]
+ x3 = p3[0]
+ y3 = p3[1]
+ z3 = p3[2]
+ x = p[0]
+ y = p[1]
+ z = p[2]
+
+ # Calculate area of triangle ABC
+ A = area(x1, y1, z1, x2, y2, z2, x3, y3, z3)
+
+ # Calculate area of triangle PBC
+ A1 = area(x, y, z, x2, y2, z2, x3, y3, z3)
+
+ # Calculate area of triangle PAC
+ A2 = area(x1, y1, z1, x, y, z, x3, y3, z3)
+
+ # Calculate area of triangle PAB
+ A3 = area(x1, y1, z1, x2, y2, z2, x, y, z)
+
+ # Check if sum of A1, A2 and A3
+ # is same as A
+ if abs(A - (A1 + A2 + A3) < 1e-6):
+ return True
+ else:
+ return False
+
+
+def intersection(planeNormal, planePoint, rayDirection, rayPoint):
+ epsilon = 1e-6
+
+ # Define plane
+ # planeNormal = np.array([0, 0, 1])
+ # planePoint = np.array([0, 0, 5]) #Any point on the plane
+
+ # Define ray
+ # rayDirection = np.array([0, -1, -1])
+ # rayPoint = np.array([0, 0, 10]) #Any point along the ray
+
+ ndotu = planeNormal.dot(rayDirection)
+
+ if abs(ndotu) < epsilon:
+ intersect = 0
+ else:
+ w = rayPoint - planePoint[0, :]
+ si = -planeNormal.dot(w) / ndotu
+ Psi = w + si * rayDirection + planePoint[0, :]
+ if isInside(planePoint[0], planePoint[1], planePoint[2], Psi) == False:
+ intersect = 0
+ else:
+ intersect = Psi[0]
+
+ return intersect
+
+intersection(np.array([0,0,1]), np.array([(1,1,0),(1,2,0),(2,1.5,0)]), np.array([0,0,1]), np.array((1.5,1.5,1)))
diff --git a/LigamentStudy/AttechmentArea.py b/LigamentStudy/AttechmentArea.py
index 653f02a..7979e8a 100644
--- a/LigamentStudy/AttechmentArea.py
+++ b/LigamentStudy/AttechmentArea.py
@@ -58,7 +58,7 @@ def cut_mesh(out2, path, i):
# select and delete vertices with negative distance
# ms4.conditional_vertex_selection(condselect="q<0.15")
- ms4.conditional_vertex_selection(condselect="q<0")
+ ms4.conditional_vertex_selection(condselect="(q <0)") # "(q <0) && (q >-0.4)"
ms4.delete_selected_vertices()
# split mesh
out4 = ms4.apply_filter('split_in_connected_components')
@@ -66,8 +66,8 @@ def cut_mesh(out2, path, i):
return ms4
-subjects = [9] # 9,13,19,23,26,29,32,35,37,41
-segments = ['tibia'] # ['femur', 'tibia'] # ['femur'] #
+subjects = [19] # 9,13,19,23,26,29,32,35,37,41
+segments = ['femur'] # ['femur', 'tibia'] # ['femur'] #
no_subjects = len(subjects)
no_segments = len(segments)
@@ -87,7 +87,6 @@ for subject in subjects:
# no_vertices = ms3.mesh(i).vertex_matrix().shape[0]
# if no_vertices < 50:
# ms3.delete_current_mesh()
- # bla
# else:
# if not os.path.isfile(path + '\Segmentation_' + segment + '_wires' + str(i-1) + '.stl') and not i == 1:
# no = i-1
@@ -107,102 +106,109 @@ for subject in subjects:
# run over all wires
error = []
mesh2 = path + 'Segmentation_' + segment + '.stl'
- ms5 = pymeshlab.MeshSet()
- ms5.load_new_mesh(mesh2)
- ms5.apply_filter('uniform_mesh_resampling', cellsize=1)
- # ms5.apply_filter('transform_rotate', rotaxis=2, angle=180)
- ms5.save_current_mesh(path + '\Segmentation_' + segment + '_resample.stl', binary=False)
-
- # for i in range(1, no_meshes+1): #range(3,4): #range(5,no_meshes+1): #
- # mesh1 = path + '\Segmentation_' + segment + '_wires' + str(i) + '.stl'
- #
- # # load meshes in new meshlab file
- # ms = pymeshlab.MeshSet()
- # ms.load_new_mesh(mesh1)
- # ms.load_new_mesh(mesh2)
- #
- # # calculate Hausdorff distance in both directions
- # out2 = ms.apply_filter('hausdorff_distance', targetmesh=1, sampledmesh=0, savesample=False, maxdist=9)
- # out1 = ms.apply_filter('hausdorff_distance', targetmesh=0, sampledmesh=1, savesample=False, maxdist=9)
- #
- # # select and delete all vertices far from the wire
- # ms.conditional_vertex_selection(condselect="q>8.9")
- # ms.delete_selected_vertices()
- # # save section containing area
- # ms.save_current_mesh(path + '\Segmentation_' + segment + '_area' + str(i) + '.stl')
- # #fit plane through section containing area
- # ms.select_all()
- # ms.fit_a_plane_to_selection()
- # plane_normal = ms.mesh(2).vertex_normal_matrix()
- #
- # ms4 = cut_mesh(out2, path, i)
- #
- # # check the number of components and remove the ones with few vertices
- # no_meshes = ms4.number_meshes()
- # meshes_to_remove = no_meshes-4
- # if meshes_to_remove > 0:
- # for ind in range(0,no_meshes):
- # no_vertices = ms4.mesh(ind).vertex_matrix().shape[0]
- # if no_vertices < 100:
- # ms4.set_current_mesh(ind)
- # ms4.delete_current_mesh()
- # else:
- # last_mesh = ind
- # else:
- # last_mesh = 3
- # # check the number of meshes
- # # if there are less than 4, split is not done in 2 large surfaces and surface needs to be closed
- # no_meshes = ms4.number_meshes()
- # if no_meshes < 4: #no_vertices < 10:
- # # load wire in new meshset
- # ms6 = pymeshlab.MeshSet()
- # ms6.load_new_mesh(path + '\Segmentation_' + segment + '_wires' + str(i) + '.stl')
- # # find for each point the largest distance on mesh
- # dist_matrix = []
- # dist_matrix_ind = []
- # start_ind = []
- # verts = ms6.mesh(0).vertex_matrix()
- # for ind in range(0, len(verts)):
- # ms6.apply_filter('colorize_by_geodesic_distance_from_a_given_point', startpoint=verts[ind], maxdistance=100)
- # dist_matrix.append(np.max(ms6.mesh(0).vertex_quality_array()))
- # dist_matrix_ind.append(np.argmax(ms6.mesh(0).vertex_quality_array()))
- # start_ind.append(ind)
- # # find which point has largest distance
- # max1 = np.argmax(dist_matrix)
- # end_point = verts[dist_matrix_ind[max1]]
- # start_point = verts[start_ind[max1]]
- # # create cylinder between these points
- # r = 0.5
- # path_cylinder = path + '\Segmentation_' + segment + '_wires' + str(i) + 'cylinder.stl'
- # cylinder_between(start_point, end_point, r, path_cylinder)
- # # combine wire and cylinder
- # ms6.load_new_mesh(path_cylinder)
- # ms6.apply_filter('mesh_boolean_union', first_mesh=0, second_mesh=1)
- # ms6.save_current_mesh(path + '\Segmentation_' + segment + '_wires' + str(i) + '.stl', binary=False)
- # # split mesh again with closed wire
- # ms4 = cut_mesh(out2, path, i)
- # # remove meshes with few vertices
- # no_meshes = ms4.number_meshes()
- # for ind in range(0,no_meshes):
- # no_vertices = ms4.mesh(ind).vertex_matrix().shape[0]
- # if no_vertices < 100:
- # ms4.set_current_mesh(ind)
- # ms4.delete_current_mesh()
- # else:
- # last_mesh = ind
- # # select last mesh to save
- # no_meshes = ms4.number_meshes()
- # # save only mesh part inside wire
- # ms4.set_current_mesh(new_curr_id=last_mesh)
- # try:
- # ms4.save_current_mesh(path + '\Segmentation_' + segment + '_area' + str(i) + '.stl', binary=False)
- # ms4.load_new_mesh(path + '\Segmentation_' + segment + '_area' + str(i) + '.stl')
- # geometric_measures = ms4.apply_filter('compute_geometric_measures')
- # surface = geometric_measures['surface_area']
- # print('Surface area ' + segment + ' ligament' + str(i) + ': ' + str(surface) + ' mm2')
- # ms4.save_project(path + '\Segmentation_' + segment + '_area' + str(i) + '.mlp')
- # except:
- # error.append(i)
+ # ms5 = pymeshlab.MeshSet()
+ # ms5.load_new_mesh(mesh2)
+ # ms5.apply_filter('uniform_mesh_resampling', cellsize=1)
+ # # ms5.apply_filter('transform_rotate', rotaxis=2, angle=180)
+ # ms5.save_current_mesh(path + '\Segmentation_' + segment + '_resample.stl', binary=False)
+
+ for i in range(3,4): #range(1, no_meshes+1): #range(3,4): #range(5,no_meshes+1): #
+ mesh1 = path + '\Segmentation_' + segment + '_wires' + str(i) + '.stl'
+
+ # load meshes in new meshlab file
+ ms = pymeshlab.MeshSet()
+ ms.load_new_mesh(mesh1)
+ ms.load_new_mesh(mesh2)
+
+ # calculate Hausdorff distance in both directions
+ out2 = ms.apply_filter('hausdorff_distance', targetmesh=1, sampledmesh=0, savesample=False, maxdist=9)
+ out1 = ms.apply_filter('hausdorff_distance', targetmesh=0, sampledmesh=1, savesample=False, maxdist=9)
+
+ # select and delete all vertices far from the wire
+ ms.conditional_vertex_selection(condselect="q>8.9")
+ ms.delete_selected_vertices()
+ # save section containing area
+ ms.save_current_mesh(path + '\Segmentation_' + segment + '_area' + str(i) + '.stl')
+ #fit plane through section containing area
+ # ms.set_current_mesh(new_curr_id=0)
+ ms.select_all()
+ ms.fit_a_plane_to_selection()
+ plane_normal = ms.mesh(2).vertex_normal_matrix()
+
+ ms4 = cut_mesh(out2, path, i)
+
+ # check the number of components and remove the ones with few vertices
+ no_meshes = ms4.number_meshes()
+ meshes_to_remove = no_meshes-4
+ if meshes_to_remove > 0:
+ for ind in range(0,no_meshes):
+ no_vertices = ms4.mesh(ind).vertex_matrix().shape[0]
+ if no_vertices < 50:
+ ms4.set_current_mesh(ind)
+ ms4.delete_current_mesh()
+ else:
+ last_mesh = ind
+ else:
+ last_mesh = 3
+ # check the number of meshes
+ # if there are less than 4, split is not done in 2 large surfaces and surface needs to be closed
+ no_meshes = ms4.number_meshes()
+ if no_meshes < 4: #no_vertices < 10:
+ # load wire in new meshset
+ ms6 = pymeshlab.MeshSet()
+ ms6.load_new_mesh(path + '\Segmentation_' + segment + '_wires' + str(i) + '.stl')
+ # find for each point the largest distance on mesh
+ dist_matrix = []
+ dist_matrix_ind = []
+ start_ind = []
+ verts = ms6.mesh(0).vertex_matrix()
+ for ind in range(0, len(verts)):
+ ms6.apply_filter('colorize_by_geodesic_distance_from_a_given_point', startpoint=verts[ind], maxdistance=100)
+ dist_matrix.append(np.max(ms6.mesh(0).vertex_quality_array()))
+ dist_matrix_ind.append(np.argmax(ms6.mesh(0).vertex_quality_array()))
+ start_ind.append(ind)
+ # find which point has largest distance
+ max1 = np.argmax(dist_matrix)
+ end_point = verts[dist_matrix_ind[max1]]
+ start_point = verts[start_ind[max1]]
+ # create cylinder between these points
+ r = 0.5
+ path_cylinder = path + '\Segmentation_' + segment + '_wires' + str(i) + 'cylinder.stl'
+ cylinder_between(start_point, end_point, r, path_cylinder)
+ # combine wire and cylinder
+ ms6.load_new_mesh(path_cylinder)
+ ms6.apply_filter('mesh_boolean_union', first_mesh=0, second_mesh=1)
+ ms6.save_current_mesh(path + '\Segmentation_' + segment + '_wires' + str(i) + '.stl', binary=False)
+ # split mesh again with closed wire
+ ms4 = cut_mesh(out2, path, i)
+ # remove meshes with few vertices
+ no_meshes = ms4.number_meshes()
+ for ind in range(0,no_meshes):
+ no_vertices = ms4.mesh(ind).vertex_matrix().shape[0]
+ if no_vertices < 50:
+ ms4.set_current_mesh(ind)
+ ms4.delete_current_mesh()
+ else:
+ last_mesh = ind
+ # select last mesh to save
+ no_meshes = ms4.number_meshes()
+ # save only mesh part inside wire
+ # ms4.set_current_mesh(new_curr_id=last_mesh)
+ to_del = [0,1,2]
+ for removes in range(len(to_del)):
+ ms4.set_current_mesh(new_curr_id=to_del[removes])
+ ms4.delete_current_mesh()
+ print(ms4.number_meshes())
+ ms4.apply_filter('flatten_visible_layers')
+ try:
+ ms4.save_current_mesh(path + '\Segmentation_' + segment + '_area' + str(i) + '.stl', binary=False)
+ ms4.load_new_mesh(path + '\Segmentation_' + segment + '_area' + str(i) + '.stl')
+ geometric_measures = ms4.apply_filter('compute_geometric_measures')
+ surface = geometric_measures['surface_area']
+ print('Surface area ' + segment + ' ligament' + str(i) + ': ' + str(surface) + ' mm2')
+ ms4.save_project(path + '\Segmentation_' + segment + '_area' + str(i) + '.mlp')
+ except:
+ error.append(i)
diff --git a/LigamentStudy/HausdorffDistance.py b/LigamentStudy/HausdorffDistance.py
new file mode 100644
index 0000000..d77679d
--- /dev/null
+++ b/LigamentStudy/HausdorffDistance.py
@@ -0,0 +1,136 @@
+import pymeshlab
+# from plyfile import PlyData, PlyElement
+import numpy as np
+import matplotlib
+import matplotlib.pyplot as plt
+# from vtk import *
+import nrrd
+import re
+import os
+import pandas as pd
+from tabulate import tabulate
+from shutil import copyfile
+import glob
+import trimesh
+
+def writeply(surface,filename):
+ """Write mesh as ply file."""
+ writer = vtkPLYWriter()
+ writer.SetInputData(surface)
+ writer.SetFileTypeToASCII()
+ writer.SetFileName(filename)
+ writer.Write()
+
+def readVTK(file):
+ reader = vtkDataSetReader()
+ reader.SetFileName(file)
+ reader.ReadAllVectorsOn()
+ reader.ReadAllScalarsOn()
+ reader.Update()
+
+ data = reader.GetOutput()
+ return data
+
+subjects = [29,32,35,37,41] #9,13,19,23,26,
+segments = ['femur'] #'femur','tibia'
+
+for segment in segments:
+ for ind, subject in enumerate(subjects):
+ 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'
+
+ # xyz_file = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone\shape_models\Segmentation_' + segment + '_' + side + '_short_' + str(
+ # subject) + reflect + '.isores.pad.com.center.aligned.clipped.cropped.tpSmoothDT_local.xyz'
+ # points1 = trimesh.load_mesh(xyz_file)
+ # # mesh = trimesh.load_mesh(path_bones + '\Segmentation_' + segment + '_' + side + '_short_' + str(subject) + '.STL')
+ # points2 = trimesh.load_mesh(r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output\femur_bone\new_bone\reconstructed\shape9.xyz')
+ # kwargs = {"scale": True}
+ # icp = trimesh.registration.icp(points2.vertices, points1.vertices, initial=np.identity(4), threshold=1e-5, max_iterations=20,**kwargs)
+ # points2.apply_transform(icp[0])
+ # # np.savetxt(r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output\femur_bone\new_bone\reconstructed\9_reconstruct_transform_icp_test.xyz', points2.vertices, delimiter=" ")
+
+ # files from SSM workflow shapeworks
+ file_com = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone\groomed\com_aligned\Segmentation_' + segment + '_' + side + '_short_' + str(
+ subject) + reflect + '.isores.pad.com.txt'
+ file_align = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone\groomed\aligned\Segmentation_' + segment + '_' + side + '_short_' + str(
+ subject) + reflect + '.isores.pad.com.center.aligned.txt'
+ pad_file = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone\groomed\padded\segementations\Segmentation_' + segment + '_' + side + '_short_' + str(
+ subject) + reflect + '.isores.pad.nrrd'
+ com_file = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone\groomed\com_aligned\Segmentation_' + segment + '_' + side + '_short_' + str(
+ subject) + reflect + '.isores.pad.com.nrrd'
+ particle_file = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone\shape_models\4096\Segmentation_' + segment + '_' + side + '_short_' + str(
+ subject) + reflect + '.isores.pad.com.center.aligned.clipped.cropped.tpSmoothDT_local.particles'
+ reconstructed_file = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone\reconstructed\mesh' + str(subject) + 'dt.stl'
+ # align_file = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output\femur_bone\new_bone\groomed\aligned\Segmentation_femur_' + side + '_short_' + str(subject) + reflect + '.isores.pad.com.center.aligned.nrrd'
+ path_bones = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone\input'
+ org_mesh = path_bones + '/Segmentation_' + segment + '_' + side + '_short_' + str(subject) + '.stl'
+
+ # get change in position from nrrd header files
+ header = nrrd.read_header(pad_file)
+ pad_position = header['space origin']
+ header = nrrd.read_header(com_file)
+ com_position = header['space origin']
+
+ # get translation from align from rotation matrix
+ rot_ssm = np.loadtxt(file_align)
+
+ # translate reconstructed SSM instance to align with original mesh
+ translate = pad_position - com_position + rot_ssm[3, :]
+ mesh3 = reconstructed_file # =local.particle file
+ ms6 = pymeshlab.MeshSet()
+ ms6.load_new_mesh(mesh3)
+ ms6.apply_filter('transform_translate_center_set_origin', traslmethod=0, axisx=translate[0],
+ axisy=translate[1], axisz=-450)
+ ms6.apply_filter('transform_translate_center_set_origin', traslmethod=0, axisx=0, axisy=0, axisz=-450)
+ ms6.apply_filter('transform_translate_center_set_origin', traslmethod=0, axisx=0, axisy=0, axisz=-450)
+ ms6.apply_filter('transform_translate_center_set_origin', traslmethod=0, axisx=0, axisy=0,
+ axisz=translate[2] + 1350)
+ ms6.save_current_mesh(path + '\SSM_' + segment + '_reconstruct_transform.stl')
+
+ # run ICP to get final position SSM point cloud on original mesh
+ mesh = trimesh.load_mesh(path + '\Segmentation_' + segment + '_resample.stl')
+ # mesh = trimesh.load_mesh(path_bones + '\Segmentation_' + segment + '_' + side + '_short_' + str(subject) + '.STL')
+ points = trimesh.load_mesh(path + '\SSM_' + segment + '_reconstruct_transform.stl')
+ if reflect == '.reflect':
+ M = trimesh.transformations.scale_and_translate((-1, 1, 1))
+ points.apply_transform(M)
+ kwargs = {"scale": False}
+ icp = trimesh.registration.icp(points.vertices, mesh, initial=np.identity(4), threshold=1e-5, max_iterations=20,**kwargs)
+ points.apply_transform(icp[0])
+ icp = trimesh.registration.icp(points.vertices, mesh, initial=np.identity(4), threshold=1e-5, max_iterations=20, **kwargs)
+ points.apply_transform(icp[0])
+ points.export(path + '\SSM_' + segment + '_reconstruct_transform_icp.stl')
+
+ ms5 = pymeshlab.MeshSet()
+ ms5.load_new_mesh(path + '\SSM_' + segment + '_reconstruct_transform_icp.stl')
+ ms5.load_new_mesh(org_mesh)
+ out1 = ms5.apply_filter('hausdorff_distance', targetmesh=1, sampledmesh=0, savesample=True)
+ out2 = ms5.apply_filter('hausdorff_distance', targetmesh=0, sampledmesh=1, savesample=True)
+
+ print('max: ' + str(max(out1['max'], out2['max'])))
+ print('min: ' + str(max(out1['min'], out2['min'])))
+ print('mean: ' + str(max(out1['mean'], out2['mean'])))
+ print('RMS: ' + str(max(out1['RMS'], out2['RMS'])))
+
+ # dist_to_use = np.argmax([out1['max'], out2['max']])
+ #
+ # vq1 = ms5.mesh(2+dist_to_use*2).vertex_quality_array()
+ #
+ # samples = [sum(vq1 < 0.5), sum((vq1 > 0.5) & (vq1 < 1)), sum((vq1 > 1) & (vq1 < 1.5)),
+ # sum((vq1 > 1.5) & (vq1 < 2)), sum(vq1 > 2)]
+ #
+ # x = np.arange(5) # the label locations
+ # width = 0.35 # the width of the bars
+ # fig, ax = plt.subplots()
+ # rects1 = ax.bar(x, samples, width, label='femoral cartilage')
+
+ ms5.save_current_mesh(path + r'/Segmentation_' + segment + '_' + side + '_short_' + str(subject) + '_HD.ply', binary=False,
+ save_vertex_quality=True)
+ np.save(path + r'/HD.np',[out1,out2])
+
+
diff --git a/LigamentStudy/OAIdownload.py b/LigamentStudy/OAIdownload.py
new file mode 100644
index 0000000..4957e53
--- /dev/null
+++ b/LigamentStudy/OAIdownload.py
@@ -0,0 +1,119 @@
+import base64
+import requests
+import json
+import urllib.request
+import shutil
+from pathlib import Path
+
+# Encode our credentials then convert it to a string.
+credentials = base64.b64encode(b'mariskawesseling:p1SM3csN5xXsGFfo').decode('utf-8')
+
+# Create the headers we will be using for all requests.
+headers = {
+ 'Authorization': 'Basic ' + credentials,
+ 'User-Agent': 'Example Client',
+ 'Accept': 'application/json'
+}
+
+# Send Http request
+response = requests.get('https://nda.nih.gov/api/package/auth', headers=headers)
+
+# Business Logic.
+
+# If the response status code does not equal 200
+# throw an exception up.
+if response.status_code != requests.codes.ok:
+ print('failed to authenticate')
+ response.raise_for_status()
+
+# The auth endpoint does no return any data to parse
+# only a Http response code is returned.
+
+# Assume code in authentication section is present.
+
+packageId = 1190875
+
+# Construct the request to get the files of package 1234
+# URL structure is: https://nda.nih.gov/api/package/{packageId}/files
+response = requests.get('https://nda.nih.gov/api/package/' + str(packageId) + '/files', headers=headers)
+
+# Get the results array from the json response.
+results = response.json()['results']
+
+# Business Logic.
+
+files = {}
+
+# Add important file data to the files dictionary.
+for f in results:
+ files[f['package_file_id']] = {'name': f['download_alias']}
+
+# Assume code in authentication section is present.
+# Assume that one of the retrieving files implementations is present too
+
+# Create a post request to the batch generate presigned urls endpoint.
+# Use keys from files dictionary to form a list, which is converted to
+# a json array which is posted.
+response = requests.post('https://nda.nih.gov/api/package/' + str(packageId) + '/files/batchGeneratePresignedUrls',
+ json=list(files.keys()), headers=headers)
+
+# Get the presigned urls from the response.
+results = response.json()['presignedUrls']
+
+# Business Logic.
+
+# Add a download key to the file's data.
+for url in results:
+ files[url['package_file_id']]['download'] = url['downloadURL']
+
+# Iterate on file id and it's data to perform the downloads.
+# for id, data in files:
+# name = data['name']
+# downloadUrl = data['download']
+# # Create a downloads directory
+# file = 'downloads/' + name
+# # Strip out the file's name for creating non-existent directories
+# directory = file[:file.rfind('/')]
+#
+# # Create non-existent directories, package files have their
+# # own directory structure, and this will ensure that it is
+# # kept in tact when downloading.
+# Path(directory).mkdir(parents=True, exist_ok=True)
+#
+# # Initiate the download.
+# with urllib.request.urlopen(downloadUrl) as dl, open(file, 'wb') as out_file:
+# shutil.copyfileobj(dl, out_file)
+
+import csv
+
+# Assume code in authentication section is present.
+
+# packageId = 1234
+
+s3Files = []
+
+# Load in and process the manifest file.
+# Not all manifest files are structured like this, all you require is
+# an S3 url and a package that has the files associated with it.
+# with open('datastructure_manifest.txt', 'r') as manifest:
+# for rows in csv.reader(manifest, dialect='excel-tab'):
+# for row in rows:
+# if row.startsWith('s3://'):
+# s3Files.append(row)
+
+# The manifest files have their column declarations listed twice, trim those out
+# s3Files = s3Files[2:]
+s3Files = ['s3://NDAR_Central_1/submission_13364/00m/0.E.1/9005075/20050926/10593811.tar.gz']
+
+# Construct the request to get the files of package 1234
+# URL structure is: https://nda.nih.gov/api/package/{packageId}/files
+response = requests.post('https://nda.nih.gov/api/package/' + str(packageId) + '/files', json=s3Files, headers=headers)
+
+# Business Logic.
+
+files = {}
+
+# Add important file data to the files dictionary.
+# We can skip having to transform the json because a json array is returned.
+for f in response.json():
+ files[f['package_file_id']] = {'name': f['download_alias']}
\ No newline at end of file
diff --git a/LigamentStudy/ProjectCentroids.py b/LigamentStudy/ProjectCentroids.py
new file mode 100644
index 0000000..4181bff
--- /dev/null
+++ b/LigamentStudy/ProjectCentroids.py
@@ -0,0 +1,213 @@
+import pandas as pd
+import os
+import trimesh
+import numpy as np
+import matplotlib.path as plt
+import copy
+import time
+
+def heron(a,b,c):
+ s = (a + b + c) / 2
+ area = (s*(s-a) * (s-b)*(s-c)) ** 0.5
+ return area
+
+def distance3d(x1,y1,z1,x2,y2,z2):
+ a=(x1-x2)**2+(y1-y2)**2 + (z1-z2)**2
+ d= a ** 0.5
+ return d
+
+def area(x1,y1,z1,x2,y2,z2,x3,y3,z3):
+ a=distance3d(x1,y1,z1,x2,y2,z2)
+ b=distance3d(x2,y2,z2,x3,y3,z3)
+ c=distance3d(x3,y3,z3,x1,y1,z1)
+ A = heron(a,b,c)
+ return A
+ # print("area of triangle is %r " %A)
+
+# A utility function to calculate area
+# of triangle formed by (x1, y1),
+# (x2, y2) and (x3, y3)
+
+# def area(x1, y1, x2, y2, x3, y3):
+# return abs((x1 * (y2 - y3) + x2 * (y3 - y1)
+# + x3 * (y1 - y2)) / 2.0)
+
+
+# A function to check whether point P(x, y)
+# lies inside the triangle formed by
+# A(x1, y1), B(x2, y2) and C(x3, y3)
+def isInside(p1, p2, p3, p):
+ x1 = p1[0]
+ y1 = p1[1]
+ z1 = p1[2]
+ x2 = p2[0]
+ y2 = p2[1]
+ z2 = p2[2]
+ x3 = p3[0]
+ y3 = p3[1]
+ z3 = p3[2]
+ x = p[0]
+ y = p[1]
+ z = p[2]
+
+ # Calculate area of triangle ABC
+ A = area(x1, y1,z1, x2, y2,z2, x3, y3,z3)
+
+ # Calculate area of triangle PBC
+ A1 = area(x, y,z, x2, y2,z2, x3, y3,z3)
+
+ # Calculate area of triangle PAC
+ A2 = area(x1, y1,z1, x, y, z,x3, y3,z3)
+
+ # Calculate area of triangle PAB
+ A3 = area(x1, y1,z1, x2, y2,z2, x, y,z)
+
+ # Check if sum of A1, A2 and A3
+ # is same as A
+ if abs(A - (A1 + A2 + A3)) < 1e-6:
+ return True
+ else:
+ return False
+
+def intersection(planeNormal,planePoint,rayDirection,rayPoint):
+ epsilon=1e-6
+
+ #Define plane
+ # planeNormal = np.array([0, 0, 1])
+ # planePoint = np.array([0, 0, 5]) #Any point on the plane
+
+ #Define ray
+ # rayDirection = np.array([0, -1, -1])
+ # rayPoint = np.array([0, 0, 10]) #Any point along the ray
+
+ ndotu = planeNormal.dot(rayDirection)
+
+ if abs(ndotu) < epsilon:
+ intersect = 0
+ else:
+ w = rayPoint - planePoint[0,:]
+ si = -planeNormal.dot(w) / ndotu
+ Psi = w + si * rayDirection + planePoint[0,:]
+ if isInside(planePoint[0], planePoint[1], planePoint[2], Psi) == False:
+ intersect = 0
+ else:
+ intersect = Psi[0]
+
+ return intersect
+
+
+subjects = [9] #[9,13,19,] #23,26,29,32,35,37,
+segment = 'femur'
+
+df = pd.read_excel(os.path.join(r"C:\Users\mariskawesseli\Documents\LigamentStudy\ImageData","surfaces2.xlsx"),
+ sheet_name='perc_of_len femur')
+lig_names = ['PCL'] #, 'MCL-p','MCL-d','posterior oblique','ACL','LCL (prox)','popliteus (dist)'
+
+for subject in subjects:
+ if subject in [9, 13, 26, 29, 32]:
+ side = 'R'
+ else:
+ side = 'L'
+ path = os.path.join(r"C:\Users\mariskawesseli\Documents\LigamentStudy\ImageData", str(subject))
+
+ mesh = path + '/Segmentation_' + segment + '_transform.stl'
+ bone = trimesh.load_mesh(mesh)
+ rot_mat = np.linalg.inv(np.loadtxt(path + '\Segmentation_' + segment + '_resample._ACS.txt'))
+
+ AP_size = np.max(bone.vertices[:,1])-np.min(bone.vertices[:,1]) # AP length
+
+ Rx = trimesh.transformations.rotation_matrix(1.57, [0, 1, 0])
+ for name in lig_names:
+ if side == 'R':
+ most_med_point = np.min(bone.vertices[:, 0]) # med
+ most_lat_point = np.max(bone.vertices[:, 0]) # lat
+ med_section_dir = + 15
+ lat_section_dir = - 15
+ else:
+ most_med_point = np.max(bone.vertices[:, 0]) # med
+ most_lat_point = np.min(bone.vertices[:, 0]) # lat
+ med_section_dir = - 10
+ lat_section_dir = + 10
+ if name == 'PCL' or name == 'ACL':
+ most_med_point = most_med_point*0.10
+ most_lat_point = most_lat_point*0.10
+ med_section = most_med_point #- med_section_dir
+ lat_section = most_lat_point #- lat_section_dir
+ else:
+ most_med_point = most_med_point
+ most_lat_point = most_lat_point
+ med_section = most_med_point + med_section_dir
+ lat_section = most_lat_point + lat_section_dir
+
+ if name == 'PCL' or name == 'MCL-p' or name == 'MCL-d' or name == 'posterior oblique':
+ LCL_points = [most_med_point * np.ones(10),
+ np.min(bone.vertices[:, 1]) + np.asarray(AP_size * df[name + 'y'][0:10]),
+ np.min(bone.vertices[:, 2]) + np.asarray(AP_size * df[name + 'z'][0:10])]
+ LCL_points = np.transpose(np.asarray(LCL_points))
+ else:
+ LCL_points = [most_lat_point*np.ones(10), np.min(bone.vertices[:,1])+np.asarray(AP_size*df[name+'y'][0:10]), np.min(bone.vertices[:,2])+np.asarray(AP_size*df[name+'z'][0:10])]
+ LCL_points = np.transpose(np.asarray(LCL_points))
+
+ for pts in range(0,10):
+ if not np.isnan(LCL_points[pts,:]).any():
+ intersect = []
+ bone_part = copy.deepcopy(bone)
+ top = max(LCL_points[:,2])+2.5
+ far_verts = bone_part.vertices[:, 2] < top
+ face_mask = far_verts[bone_part.faces].all(axis=1)
+ bone_part.update_faces(face_mask)
+ if name == 'PCL' or name == 'MCL-p' or name == 'MCL-d' or name == 'posterior oblique':
+ if side == 'R':
+ far_verts = bone_part.vertices[:, 0] < med_section
+ else:
+ far_verts = bone_part.vertices[:, 0] > med_section
+ face_mask = far_verts[bone_part.faces].all(axis=1)
+ bone_part.update_faces(face_mask)
+ # trimesh.Scene(bone_part).show()
+ else:
+ if side == 'R':
+ far_verts = bone_part.vertices[:, 0] > lat_section
+ else:
+ far_verts = bone_part.vertices[:, 0] < lat_section
+ face_mask = far_verts[bone_part.faces].all(axis=1)
+ bone_part.update_faces(face_mask)
+ # trimesh.Scene(bone_part).show()
+ # tic = time.perf_counter()
+ for count, tr in enumerate(bone_part.face_normals):
+ intersect.append(intersection(tr, bone_part.vertices[bone_part.faces[count,:]], np.array([1,0,0]), LCL_points[pts,:]))
+ # toc = time.perf_counter()
+ # print(f"Downloaded the tutorial in {toc - tic:0.4f} seconds")
+
+ # T = trimesh.transformations.translation_matrix(LCL_points[pts])
+ # point = trimesh.creation.cylinder(0.5, height=0.5, sections=None, segment=None, transform=T)
+ # trimesh.Scene([bone_part, point]).show()
+
+ x_coord = [i for i in intersect if i != 0]
+ if not len(x_coord) == 0:
+ if name == 'MCL-p' or name == 'MCL-d' or name == 'posterior oblique':
+ to_use = np.argmin(abs(x_coord - most_med_point))
+ elif name == 'PCL':
+ to_use = np.argmin(abs(x_coord - most_med_point))
+ elif name == 'ACL':
+ to_use = np.argmin(abs(x_coord - most_lat_point))
+ else:
+ to_use = np.argmax(abs(x_coord - most_lat_point))
+ if not abs(x_coord[to_use]-LCL_points[pts, 0]) > 20:
+ LCL_points[pts, 0] = x_coord[to_use]
+
+
+ # points = trimesh.PointCloud(LCL_points, colors=None, metadata=None) # create point cloud
+
+ points = []
+ for ind in range(0,10):
+ T = trimesh.transformations.translation_matrix(LCL_points[ind])
+ R = np.linalg.inv(rot_mat)
+ M = trimesh.transformations.concatenate_matrices(R, T, Rx)
+ point = trimesh.creation.cylinder(0.5, height=0.5, sections=None, segment=None, transform=M)
+ # point = trimesh.creation.icosphere(subdivisions=3, radius=1.0, color=None, transform=T)
+ if ind == 0:
+ points = point
+ else:
+ points = trimesh.boolean.union([points,point])
+
+ points.export(os.path.join(r"C:\Users\mariskawesseli\Documents\LigamentStudy\ImageData", str(subject),name+'centroids.stl'))
diff --git a/LigamentStudy/SlicerExportXray.py b/LigamentStudy/SlicerExportXray.py
new file mode 100644
index 0000000..4f2af35
--- /dev/null
+++ b/LigamentStudy/SlicerExportXray.py
@@ -0,0 +1,173 @@
+import glob
+import shutil
+import os
+import DICOMScalarVolumePlugin
+import slicer
+import vtk
+#exec(open(r'C:\Users\mariskawesseli\Documents\GitLab\Other\LigamentStudy\SlicerExportXray.py').read())
+
+subjects = [13,19,23,26,29,32,35,37,41] # 9
+for subject in subjects:
+ lig_names = ['PCL', 'MCL-p','MCL-d','posterior oblique','ACL','LCL (prox)','popliteus (dist)']
+ path = os.path.join(r"C:\Users\mariskawesseli\Documents\LigamentStudy\ImageData", str(subject),'DRR')
+ slicer.mrmlScene.Clear(0)
+ slicer.util.loadScene(glob.glob(os.path.join( path,"*.mrml"))[0])
+ no_med=-1
+ no_lat=-1
+ for name in lig_names:
+ slicer.util.loadSegmentation(os.path.join(r'C:\Users\mariskawesseli\Documents\LigamentStudy\ImageData',str(subject),name+'centroids.stl'))
+
+ if name == 'PCL' or name == 'MCL-p' or name == 'MCL-d' or name == 'posterior oblique':
+ segmentationNode = slicer.util.getNode('Segmentation_med')
+ else:
+ segmentationNode = slicer.util.getNode('Segmentation_lat')
+ segmentationNode.GetSegmentation().CopySegmentFromSegmentation(slicer.util.getNode(name+'centroids').GetSegmentation(),name+'centroids')
+
+ labelmapVolumeNode = slicer.mrmlScene.AddNewNodeByClass('vtkMRMLLabelMapVolumeNode')
+ # slicer.modules.segmentations.logic().ExportVisibleSegmentsToLabelmapNode(segmentationNode, labelmapVolumeNode)
+ segmentIds = vtk.vtkStringArray()
+ segmentIds.InsertNextValue(name + 'centroids')
+ slicer.vtkSlicerSegmentationsModuleLogic.ExportSegmentsToLabelmapNode(segmentationNode, segmentIds, labelmapVolumeNode)
+
+ outputvolumenode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLScalarVolumeNode", 'Labelmap'+name)
+ sef = slicer.modules.volumes.logic().CreateScalarVolumeFromVolume(slicer.mrmlScene, outputvolumenode, labelmapVolumeNode)
+ volumeNode = slicer.util.getNode("Labelmap"+name)
+ voxels = slicer.util.arrayFromVolume(volumeNode)
+ voxels[voxels==1] = 8000
+ voxels[voxels==2] = 8000
+ voxels[voxels==3] = 8000
+ voxels[voxels==4] = 8000
+ voxels[voxels==0] = -8000
+
+ rtImagePlan = slicer.util.getNode("RTPlan")
+ if name=='PCL' or name=='MCL-p' or name == 'MCL-d' or name=='posterior oblique':
+ beam_name = "NewBeam_med"
+ no_med +=1
+ no=no_med
+ else:
+ beam_name = "NewBeam_lat"
+ no_lat +=1
+ no=no_lat
+ rtImageBeam = rtImagePlan.GetBeamByName(beam_name)
+ Volume = slicer.util.getNode("Labelmap"+name)
+ # Create DRR image computation node for user imager parameters
+ drrParameters = slicer.mrmlScene.AddNewNodeByClass('vtkMRMLDrrImageComputationNode', 'rtImageBeamParams')
+ # Set and observe RTImage beam by the DRR node
+ drrParameters.SetAndObserveBeamNode(rtImageBeam)
+ # Get DRR computation logic
+ drrLogic = slicer.modules.drrimagecomputation.logic()
+ # Update imager markups for the 3D view and slice views (optional)
+ drrLogic.UpdateMarkupsNodes(drrParameters)
+ # Update imager normal and view-up vectors (mandatory)
+ drrLogic.UpdateNormalAndVupVectors(drrParameters) # REQUIRED
+ # Compute DRR image
+ drr_image = drrLogic.ComputePlastimatchDRR(drrParameters, Volume)
+ # slicer.mrmlScene.Clear(0)
+ if no == 0:
+ volumeNode = slicer.util.getNode("DRR : " + beam_name)
+ else:
+ volumeNode = slicer.util.getNode("DRR : " + beam_name + "_" + str(no))
+ outputFolder = os.path.join(r"C:\Users\mariskawesseli\Documents\LigamentStudy\ImageData", str(subject), "DRR")
+ # Create patient and study and put the volume under the study
+ shNode = slicer.vtkMRMLSubjectHierarchyNode.GetSubjectHierarchyNode(slicer.mrmlScene)
+ patientItemID = shNode.CreateSubjectItem(shNode.GetSceneItemID(), "test patient")
+ studyItemID = shNode.CreateStudyItem(patientItemID, "test study")
+ volumeShItemID = shNode.GetItemByDataNode(volumeNode)
+ shNode.SetItemParent(volumeShItemID, studyItemID)
+ exporter = DICOMScalarVolumePlugin.DICOMScalarVolumePluginClass()
+ exportables = exporter.examineForExport(volumeShItemID)
+ for exp in exportables:
+ exp.directory = outputFolder
+
+ exporter.export(exportables)
+ folders = [x[0] for x in os.walk(outputFolder)]
+ im_folder = [s for s in folders if str(volumeShItemID) in s]
+ shutil.move(im_folder[0] + '\IMG0001.dcm', outputFolder + '/' + name + '0001.dcm')
+ os.rmdir(im_folder[0])
+
+ names = ['all','med','lat']
+ for name in names:
+ volumeNode = slicer.util.getNode("Segmentation_"+name+'-label')
+ voxels = slicer.util.arrayFromVolume(volumeNode)
+ voxels[voxels==1] = 8000
+ voxels[voxels==2] = 8000
+ voxels[voxels==3] = 8000
+ voxels[voxels==4] = 8000
+ voxels[voxels == 5] = 8000
+ voxels[voxels == 6] = 8000
+ voxels[voxels == 7] = 8000
+ voxels[voxels == 8] = 8000
+ voxels[voxels==0] = -8000
+
+ names = ["med_fem", "lat_fem", "med_wires", "lat_wires", "med_all_wires", "lat_all_wires"]
+ for name in names:
+ rtImagePlan = slicer.util.getNode("RTPlan")
+ if 'lat' in name:
+ beam_name = "NewBeam_lat"
+ no_lat += 1
+ no = no_lat
+ else:
+ beam_name = "NewBeam_med"
+ no_med += 1
+ no = no_med
+ rtImageBeam = rtImagePlan.GetBeamByName(beam_name)
+ if 'fem' in name:
+ Volume = slicer.util.getNode("resampled06")
+ elif 'med_wires' in name:
+ Volume = slicer.util.getNode("Segmentation_med-label")
+ elif 'lat_wires' in name:
+ Volume = slicer.util.getNode("Segmentation_lat-label")
+ elif 'all' in name:
+ Volume = slicer.util.getNode("Segmentation_all-label")
+ # Create DRR image computation node for user imager parameters
+ drrParameters = slicer.mrmlScene.AddNewNodeByClass('vtkMRMLDrrImageComputationNode', 'rtImageBeamParams')
+ # Set and observe RTImage beam by the DRR node
+ drrParameters.SetAndObserveBeamNode(rtImageBeam)
+ # Get DRR computation logic
+ drrLogic = slicer.modules.drrimagecomputation.logic()
+ # Update imager markups for the 3D view and slice views (optional)
+ drrLogic.UpdateMarkupsNodes(drrParameters)
+ # Update imager normal and view-up vectors (mandatory)
+ drrLogic.UpdateNormalAndVupVectors(drrParameters) # REQUIRED
+ # Compute DRR image
+ drr_image = drrLogic.ComputePlastimatchDRR(drrParameters, Volume)
+ # slicer.mrmlScene.Clear(0)
+
+ if no == 0:
+ volumeNode = slicer.util.getNode("DRR : " + beam_name) #getNode("DRR : Beam_" + name)
+ else:
+ volumeNode = slicer.util.getNode("DRR : " + beam_name + '_' + str(no))
+ outputFolder = os.path.join(r"C:\Users\mariskawesseli\Documents\LigamentStudy\ImageData", str(subject), "DRR")
+ # Create patient and study and put the volume under the study
+ shNode = slicer.vtkMRMLSubjectHierarchyNode.GetSubjectHierarchyNode(slicer.mrmlScene)
+ patientItemID = shNode.CreateSubjectItem(shNode.GetSceneItemID(), "test patient")
+ studyItemID = shNode.CreateStudyItem(patientItemID, "test study")
+ volumeShItemID = shNode.GetItemByDataNode(volumeNode)
+ shNode.SetItemParent(volumeShItemID, studyItemID)
+ exporter = DICOMScalarVolumePlugin.DICOMScalarVolumePluginClass()
+ exportables = exporter.examineForExport(volumeShItemID)
+ for exp in exportables:
+ exp.directory = outputFolder
+
+ exporter.export(exportables)
+ folders = [x[0] for x in os.walk(outputFolder)]
+ im_folder = [s for s in folders if str(volumeShItemID) in s]
+ shutil.move(im_folder[0] + '\IMG0001.dcm', outputFolder+'/' + name + '0001.dcm')
+ os.rmdir(im_folder[0])
+
+# in slicer
+# import resampled data
+# import segmented seperate wires as segmentation
+# create 3 new segmentations (all, med, lat) with resampled image as master volume
+# in segmentations - copy wires segmentation to segmentation resampled volume
+# add correct wires to med/lat/all
+# export visible segments to label map
+# in volumes - convert to scalar volume
+
+# import resampled femur
+# external beam planning
+# Ref volume: resampled06
+# Gantry: 101/281
+# Structure set: segmentsation all
+# DRR image computation
+# export to DICOM - crate dicom series
diff --git a/LigamentStudy/SlicerPositionBeam.py b/LigamentStudy/SlicerPositionBeam.py
new file mode 100644
index 0000000..5a903be
--- /dev/null
+++ b/LigamentStudy/SlicerPositionBeam.py
@@ -0,0 +1,36 @@
+#exec(open(r'C:\Users\mariskawesseli\Documents\GitLab\Other\LigamentStudy\SlicerPositionBeam.py').read())
+import os,glob
+subject = 23
+
+path = os.path.join(r"C:\Users\mariskawesseli\Documents\LigamentStudy\ImageData", str(subject),'DRR')
+slicer.mrmlScene.Clear(0)
+slicer.util.loadScene(glob.glob(os.path.join(path,"*.mrml"))[0])
+# Create dummy RTPlan
+rtImagePlan = getNode("RTPlan")
+rtImageBeam = rtImagePlan.GetBeamByName("NewBeam_lat")
+# Set required beam parameters
+current_angle = rtImageBeam.GetGantryAngle()
+rtImageBeam.SetGantryAngle(current_angle-7)
+rtImageBeam.SetCouchAngle(355)
+
+rtImageBeam2 = rtImagePlan.GetBeamByName("NewBeam_med")
+rtImageBeam2.SetGantryAngle(current_angle-7+180)
+rtImageBeam2.SetCouchAngle(355)
+
+# # Get CT volume
+# ctVolume = getNode('resampled06')
+# # Create DRR image computation node for user imager parameters
+# drrParameters = slicer.mrmlScene.AddNewNodeByClass('vtkMRMLDrrImageComputationNode', 'rtImageBeamParams')
+# # Set and observe RTImage beam by the DRR node
+# drrParameters.SetAndObserveBeamNode(rtImageBeam)
+# # Get DRR computation logic
+# drrLogic = slicer.modules.drrimagecomputation.logic()
+# # Update imager markups for the 3D view and slice views (optional)
+# drrLogic.UpdateMarkupsNodes(drrParameters)
+# # Update imager normal and view-up vectors (mandatory)
+# drrLogic.UpdateNormalAndVupVectors(drrParameters) # REQUIRED
+# # Compute DRR image
+# drrLogic.ComputePlastimatchDRR(drrParameters, ctVolume)
+
+# save scene
+slicer.util.saveScene(glob.glob(os.path.join(path,"*.mrml"))[0])
diff --git a/LigamentStudy/VisualiseSSM.py b/LigamentStudy/VisualiseSSM.py
index abcea63..b8f1f96 100644
--- a/LigamentStudy/VisualiseSSM.py
+++ b/LigamentStudy/VisualiseSSM.py
@@ -1,9 +1,10 @@
-import vtk
import sys
import os
import vtk
from numpy import random, genfromtxt, size
import trimesh
+import numpy as np
+
class VtkPointCloud:
def __init__(self, zMin=-10.0, zMax=10.0, maxNumPoints=1e6):
@@ -41,6 +42,7 @@ class VtkPointCloud:
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)):
@@ -49,6 +51,7 @@ def load_data(data, pointCloud):
return pointCloud
+
def load_stl(filename):
reader = vtk.vtkSTLReader()
reader.SetFileName(filename)
@@ -64,6 +67,7 @@ def load_stl(filename):
return actor
+
def create_pointcloud_polydata(points, colors=None):
"""https://github.com/lmb-freiburg/demon
Creates a vtkPolyData object with the point cloud from numpy arrays
@@ -89,7 +93,8 @@ def create_pointcloud_polydata(points, colors=None):
vcolors.SetName("Colors")
vcolors.SetNumberOfTuples(points.shape[0])
for i in range(points.shape[0]):
- vcolors.SetTuple3(i, colors[0], colors[1], colors[2])
+ c = rgb(1,10,colors[i])
+ vcolors.SetTuple3(i, c[0], c[1], c[2])
vpoly.GetPointData().SetScalars(vcolors)
vcells = vtk.vtkCellArray()
@@ -103,35 +108,53 @@ def create_pointcloud_polydata(points, colors=None):
return vpoly
+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
+
+
if __name__ == '__main__':
- subjects = [13] #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],
- [6,5,6,6,6,6,4,4,5,5],
- [3,2,5,3,3,2,2,0,3,3],
- [7,8,7,7,7,5,7,6,7,0],
- [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],
- [0,3,8,0,0,0,0,0,0,0]]
- ligaments_tib = [[5,7,6,5,3,4,4,5,5,4],
- [3,3,7,3,5,3,5,4,3,3],
- [1,1,1,1,1,1,1,1,1,1],
- [4,5,3,4,4,5,3,2,4,3],
- [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],
- [0,0,0,0,0,0,0,0,0,0]]
+ subjects = ['S0'] #9,13,19,23,26,29,32,35,37,41
+
+ 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],
+ [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]]
for segment in segments:
- SSMpoints = [[] for i in range(8)]
- for ind in range(0,8):
+ SSMpoints = [[] for i in range(11)]
+ 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')
+ path = os.path.join(r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone_short\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))
@@ -150,24 +173,36 @@ if __name__ == '__main__':
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 + '\meanshape_ligs_color.xyz')[:, 3]
+ # point_cloud_lig = create_pointcloud_polydata(points_lig, color)
bone_actor = load_stl(path + '/mean_shape.stl')
bone_actor.GetProperty().SetOpacity(0.75)
else:
- points_lig = trimesh.load_mesh(path + '\SSM_' + segment + '_pred_points.xyz') #_areas
+ points_lig = trimesh.load_mesh(path + '\SSM_' + segment + '_short_pred_points.xyz') #_pred_points_color
point_cloud_lig = create_pointcloud_polydata(points_lig)
- bone_actor = load_stl(path + '/Segmentation_' + segment + '_resample.stl')
+ # 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 _pred_points
+ # point_cloud_lig = create_pointcloud_polydata(points_lig,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')
+ else:
+ bone_actor = load_stl(path + '/Segmentation_' + segment + '_resample.stl')
+ wire_actor = load_stl(path + '/Segmentation_' + segment + '_wires.stl')
+ wire_actor.GetProperty().SetColor(1, 1, 0)
bone_actor.GetProperty().SetOpacity(0.75)
- wire_actor = load_stl(path + '/Segmentation_' + segment + '_wires.stl')
- wire_actor.GetProperty().SetColor(0, 0, 1)
points_bone = trimesh.load_mesh(path + '\SSM_' + segment + '_transform_icp.xyz')
point_cloud_bone = create_pointcloud_polydata(points_bone)
- mapper = vtk.vtkPolyDataMapper()
- mapper.SetInputData(point_cloud_bone)
- actor = vtk.vtkActor()
- actor.SetMapper(mapper)
- actor.GetProperty().SetColor(0,0,0)
- actor.GetProperty().SetPointSize(2)
+ # 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)
mapper2 = vtk.vtkPolyDataMapper()
mapper2.SetInputData(point_cloud_lig)
@@ -181,7 +216,7 @@ if __name__ == '__main__':
# renderer.AddActor(actor)
renderer.AddActor(actor2)
renderer.AddActor(bone_actor)
- if not subject==100:
+ if not subject == 100 and not subject == 'S0':
renderer.AddActor(wire_actor)
# renderer.SetBackground(.2, .3, .4)
renderer.SetBackground(1.0, 1.0, 1.0)
diff --git a/LigamentStudy/VisualizeProjectedCentroids.py b/LigamentStudy/VisualizeProjectedCentroids.py
new file mode 100644
index 0000000..99d8876
--- /dev/null
+++ b/LigamentStudy/VisualizeProjectedCentroids.py
@@ -0,0 +1,207 @@
+import vtk
+import sys
+import os
+import vtk
+from numpy import random, genfromtxt, size
+import trimesh
+
+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):
+ """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])
+ vpoly = vtk.vtkPolyData()
+ vpoly.SetPoints(vpoints)
+
+ if not colors is None:
+ vcolors = vtk.vtkUnsignedCharArray()
+ vcolors.SetNumberOfComponents(3)
+ vcolors.SetName("Colors")
+ vcolors.SetNumberOfTuples(points.shape[0])
+ for i in range(points.shape[0]):
+ vcolors.SetTuple3(i, colors[0], colors[1], colors[2])
+ vpoly.GetPointData().SetScalars(vcolors)
+
+ vcells = vtk.vtkCellArray()
+
+ for i in range(points.shape[0]):
+ vcells.InsertNextCell(1)
+ vcells.InsertCellPoint(i)
+
+ vpoly.SetVerts(vcells)
+
+ return vpoly
+
+
+lig_names = ['PCL', 'MCL-p','MCL-d','posterior oblique','ACL','LCL (prox)','popliteus (dist)']
+color = ((0.75,1,0.5),
+ (0,0.5,0),
+ (1,0,1),
+ (0.5,1,1),
+ (1,1,0),
+ (1,0.5,0.75),
+ (1,0.5,0))
+
+if __name__ == '__main__':
+ subjects = [9] #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],
+ [6,5,6,6,6,6,4,4,5,5],
+ [3,2,5,3,3,2,2,0,3,3],
+ [7,8,7,7,7,5,7,6,7,0],
+ [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],
+ [0,3,8,0,0,0,0,0,0,0]]
+ ligaments_tib = [[5,7,6,5,3,4,4,5,5,4],
+ [3,3,7,3,5,3,5,4,3,3],
+ [1,1,1,1,1,1,1,1,1,1],
+ [4,5,3,4,4,5,3,2,4,3],
+ [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],
+ [0,0,0,0,0,0,0,0,0,0]]
+
+ for segment in segments:
+ SSMpoints = [[] for i in range(8)]
+ for ind in range(0,8):
+ 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')
+ 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)
+ bone_actor = load_stl(path + '/mean_shape.stl')
+ bone_actor.GetProperty().SetOpacity(0.75)
+ else:
+ points_lig = trimesh.load_mesh(path + '\SSM_' + segment + '_pred_points.xyz') # _areas
+ point_cloud_lig = create_pointcloud_polydata(points_lig)
+ bone_actor = load_stl(path + '/Segmentation_' + segment + '_resample.stl')
+ bone_actor.GetProperty().SetOpacity(0.75)
+ wire_actor = load_stl(path + '/Segmentation_' + segment + '_wires.stl')
+ wire_actor.GetProperty().SetColor(0, 0, 1)
+ lig_actor = []
+ for count, lig in enumerate(lig_names):
+ lig_actor.append(load_stl(os.path.join(path,lig+'centroids.stl')))
+ lig_actor[count].GetProperty().SetColor(color[count])
+
+ mapper2 = vtk.vtkPolyDataMapper()
+ mapper2.SetInputData(point_cloud_lig)
+ actor2 = vtk.vtkActor()
+ actor2.SetMapper(mapper2)
+ actor2.GetProperty().SetColor(1, 0, 0)
+ actor2.GetProperty().SetPointSize(5)
+
+ # Renderer
+ renderer = vtk.vtkRenderer()
+ renderer.AddActor(bone_actor)
+ if not subject==100:
+ renderer.AddActor(wire_actor)
+ for count, lig in enumerate(lig_names):
+ renderer.AddActor(lig_actor[count])
+ # renderer.AddActor(actor2)
+ renderer.SetBackground(1.0, 1.0, 1.0)
+ renderer.ResetCamera()
+
+ # Render Window
+ renderWindow = vtk.vtkRenderWindow()
+ renderWindow.AddRenderer(renderer)
+
+ # Interactor
+ renderWindowInteractor = vtk.vtkRenderWindowInteractor()
+ renderWindowInteractor.SetRenderWindow(renderWindow)
+ renderWindowInteractor.GetInteractorStyle().SetCurrentStyleToTrackballCamera()
+
+ # Begin Interaction
+ renderWindow.Render()
+ renderWindow.SetWindowName("XYZ Data Viewer")
+ renderWindowInteractor.Start()
+
diff --git a/LigamentStudy/Xray.py b/LigamentStudy/Xray.py
index 97645d0..182c281 100644
--- a/LigamentStudy/Xray.py
+++ b/LigamentStudy/Xray.py
@@ -1,71 +1,101 @@
-# in slicer
-# import resampled data
-# import segmented seperate wires as segmentation
-# create 3 new segmentations (all, med, lat) with resampled image as master volume
-# in segmentations - copy wires segmentation to segmentation resampled volume
-# add correct wires to med/lat/all
-# export visible segments to label map
-# in volumes - convert to scalar volume
-# volumeNode = getNode("Segmentation_lat-label")
-# voxels = slicer.util.arrayFromVolume(volumeNode)
-# voxels[voxels==1] = 8000
-# voxels[voxels==0] = -8000
-# import resampled femur
-# external beam planning
-# Ref volume: resampled06
-# Gantry: 101/281
-# Structure set: segmentsation all
-# DRR image computation
-# export to DICOM - crate dicom series
-
import cv2
import numpy as np
import SimpleITK as sitk
import pydicom as dicom
import os
+import glob
+
+subjects = [9,13,19,23,26,29,32,35,37,41]
+for subject in subjects:
+ path_drr = r'C:/Users/mariskawesseli/Documents/LigamentStudy/ImageData/'+str(subject)+'/DRR/'
+ images = ['med_fem0001.dcm','med_wires0001.dcm','lat_fem0001.dcm','lat_wires0001.dcm',
+ 'med_fem0001.dcm','med_all_wires0001.dcm','lat_fem0001.dcm','lat_all_wires0001.dcm']
+ lig_names = ['PCL', 'MCL-p','MCL-d','posterior oblique','ACL','LCL (prox)','popliteus (dist)']
+
+ for file in glob.glob(os.path.join(path_drr,"*.dcm")):
+ image = file
+ ds = dicom.dcmread(image)
+ pixel_array_numpy = ds.pixel_array
+ image = image.replace('.dcm', '.jpg')
+ cv2.imwrite(os.path.join(path_drr, image), pixel_array_numpy)
-path_drr = r'C:/Users/mariskawesseli/Documents/LigamentStudy/ImageData/41/DRR/'
-images = ['med_fem0001.dcm','med_wires0001.dcm','lat_fem0001.dcm','lat_wires0001.dcm',
- 'med_fem0001.dcm','med_all_wires0001.dcm','lat_fem0001.dcm','lat_all_wires0001.dcm']
-for im in images:
- image = path_drr + '/' + im
- ds = dicom.dcmread(image)
- pixel_array_numpy = ds.pixel_array
- image = image.replace('.dcm', '.jpg')
- cv2.imwrite(os.path.join(path_drr, image), pixel_array_numpy)
+ for ind in range(0,4):
+ if subject in [9, 13, 26, 29, 32]:
+ side = 'R'
+ else:
+ side = 'L'
-for ind in range(0,4):
- src = cv2.imread(os.path.join(path_drr,images[0+2*ind].replace('.dcm', '.jpg')))
- mask = cv2.imread(os.path.join(path_drr, images[1+2*ind].replace('.dcm', '.jpg')))
+ mask2 = []
+ src = cv2.imread(os.path.join(path_drr,images[0+2*ind].replace('.dcm', '.jpg')))
+ mask = cv2.imread(os.path.join(path_drr, images[1+2*ind].replace('.dcm', '.jpg')))
+ if 'med' in images[0+2*ind]:
+ for ind2 in range(0,4):
+ mask2.append(cv2.imread(os.path.join(path_drr, lig_names[ind2] + '0001.dcm'.replace('.dcm', '.jpg'))))
+ else:
+ for ind2 in range(4, 7):
+ mask2.append(cv2.imread(os.path.join(path_drr, lig_names[ind2] + '0001.dcm'.replace('.dcm', '.jpg'))))
- # convert mask to gray and then threshold it to convert it to binary
- gray = cv2.cvtColor(mask, cv2.COLOR_BGR2GRAY)
- # ret, binary = cv2.threshold(gray, 40, 255, cv2.THRESH_BINARY)
- blur = cv2.GaussianBlur(gray, (3,3), 0)
- binary = cv2.threshold(blur, 250, 255, cv2.THRESH_BINARY_INV)[1] # + cv2.THRESH_OTSU
+ # convert mask to gray and then threshold it to convert it to binary
+ gray = cv2.cvtColor(mask, cv2.COLOR_BGR2GRAY)
+ gray2 = []
+ for new_mask in mask2:
+ gray2.append(cv2.cvtColor(new_mask, cv2.COLOR_BGR2GRAY))
+ # ret, binary = cv2.threshold(gray, 40, 255, cv2.THRESH_BINARY)
+ blur = cv2.GaussianBlur(gray, (3,3), 0)
+ blur2 = []
+ for new_gray in gray2:
+ blur2.append(cv2.GaussianBlur(new_gray, (3, 3), 0))
+ binary = cv2.threshold(blur, 250, 255, cv2.THRESH_BINARY_INV)[1] # + cv2.THRESH_OTSU
+ binary2 = []
+ for new_blur in blur2:
+ binary2.append(cv2.threshold(new_blur, 250, 255, cv2.THRESH_BINARY_INV)[1]) # + cv2.THRESH_OTSU
- # find contours of two major blobs present in the mask
- contours,hierarchy = cv2.findContours(binary, cv2.RETR_LIST, cv2.CHAIN_APPROX_NONE)
+ # find contours of two major blobs present in the mask
+ contours,hierarchy = cv2.findContours(binary, cv2.RETR_LIST, cv2.CHAIN_APPROX_NONE)
+ contours2 = []
+ for new_binary in binary2:
+ contoursX, hierarchyX = cv2.findContours(new_binary, cv2.RETR_LIST, cv2.CHAIN_APPROX_NONE)
+ contours2.append(contoursX)
- # draw the found contours on to source image
- for contour in contours:
- cv2.drawContours(src, contour, -1, (255,0,0), thickness = 1)
+ # draw the found contours on to source image
+ for contour in contours:
+ cv2.drawContours(src, contour, -1, (255,0,0), thickness = 1)
+ colors = [(0,0,255),(0,255,0),(255,0,255),(0,255,255)]
+ count = -1
+ for new_contours in contours2:
+ count += 1
+ c = colors[count]
+ for contour in new_contours:
+ cv2.drawContours(src, contour, -1, c, thickness = 1)
- # split source to B,G,R channels
- b,g,r = cv2.split(src)
+ # split source to B,G,R channels
+ b,g,r = cv2.split(src)
+ b2, g2, r2 = cv2.split(src)
- # add a constant to R channel to highlight the selected area in red
- r = cv2.add(b, 30, dst = b, mask = binary, dtype = cv2.CV_8U)
+ # add a constant to R channel to highlight the selected area in red
+ r = cv2.add(b, 30, dst = b, mask = binary, dtype = cv2.CV_8U)
+ #r2 = []
+ # for new_binary in binary2:
+ # r2 = cv2.add(b, 30, dst=b, mask=new_binary, dtype=cv2.CV_8U)
- # merge the channels back together
- img_overlay = cv2.merge((b,g,r), src)
- # cv2.imshow('overlay', img_overlay)
- if (ind == 0) or (ind == 2):
- img_rot = cv2.rotate(img_overlay, cv2.ROTATE_180)
- img_rot = cv2.flip(img_rot, 1)
- else:
- img_rot = img_overlay
- cv2.imwrite(os.path.join(path_drr, images[1+2*ind].replace('.dcm', '_combine.jpg')),img_rot)
+ # merge the channels back together
+ img_overlay = cv2.merge((b,g,r), src)
+ # for new_r in r2:
+ img_overlay = cv2.merge((b2, g2, r2), img_overlay)
+ # cv2.imshow('overlay', img_overlay)
+ if side == 'R':
+ if (ind == 1) or (ind == 3):
+ img_rot = cv2.rotate(img_overlay, cv2.ROTATE_180)
+ img_rot = cv2.flip(img_rot, 1)
+ else:
+ img_rot = img_overlay
+ else:
+ if (ind == 0) or (ind == 2):
+ img_rot = cv2.rotate(img_overlay, cv2.ROTATE_180)
+ img_rot = cv2.flip(img_rot, 1)
+ else:
+ img_rot = img_overlay
+ cv2.imwrite(os.path.join(path_drr, images[1+2*ind].replace('.dcm', '_combine.jpg')),img_rot)
diff --git a/LigamentStudy/fitSSM.py b/LigamentStudy/fitSSM.py
index 9e674e8..1cc791b 100644
--- a/LigamentStudy/fitSSM.py
+++ b/LigamentStudy/fitSSM.py
@@ -9,24 +9,32 @@ from tabulate import tabulate
from shutil import copyfile
import glob
-subjects = [9,13,19,23,26,29,32,35,37,41]
-segments = ['femur'] #'femur',
+subjects = [9,13,19,23,26,29,32,35,37,41] #
+segments = ['tibia'] #'femur',
+short = 1
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],
- [7,8,7,7,7,5,7,6,7,0],
- [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],
- [0,3,8,0,0,0,0,0,0,0]]
+ [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],
- [3,3,7,3,5,3,5,4,3,3],
- [1,1,1,1,1,1,1,1,1,1],
- [4,5,3,4,4,5,3,2,4,3],
- [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],
- [0,0,0,0,0,0,0,0,0,0]]
+ [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]]
for segment in segments:
if segment == 'femur':
@@ -34,8 +42,8 @@ for segment in segments:
else:
ligaments = ligaments_tib
- SSMpoints = [[] for i in range(8)]
- for ind in range(0,8):
+ SSMpoints = [[] for i in range(11)]
+ for ind in range(0,11):
SSMpoints[ind] = [[] for i in range(10)]
for ind, subject in enumerate(subjects):
@@ -49,69 +57,86 @@ for segment in segments:
"""SSM part"""
# files from SSM workflow shapeworks
- file_com = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone\groomed\com_aligned\Segmentation_' + segment + '_' + side + '_short_' + str(subject) + reflect + '.isores.pad.com.txt'
- file_align = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone\groomed\aligned\Segmentation_' + segment + '_' + side + '_short_' + str(subject) + reflect + '.isores.pad.com.center.aligned.txt'
- pad_file = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone\groomed\padded\segementations\Segmentation_' + segment + '_' + side + '_short_' + str(subject) + reflect + '.isores.pad.nrrd'
- com_file = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone\groomed\com_aligned\Segmentation_' + segment + '_' + side + '_short_' + str(subject) + reflect + '.isores.pad.com.nrrd'
- particle_file = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone\shape_models\4096\Segmentation_' + segment + '_' + side + '_short_' + str(subject) + reflect + '.isores.pad.com.center.aligned.clipped.cropped.tpSmoothDT_local.particles'
- xyz_file = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone\shape_models\Segmentation_' + segment + '_' + side + '_short_' + str(subject) + reflect + '.isores.pad.com.center.aligned.clipped.cropped.tpSmoothDT_local.xyz'
- # align_file = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output\femur_bone\new_bone\groomed\aligned\Segmentation_femur_' + side + '_short_' + str(subject) + reflect + '.isores.pad.com.center.aligned.nrrd'
- path_bones = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone\input'
-
- # # get change in position from nrrd header files
- # header = nrrd.read_header(pad_file)
- # pad_position = header['space origin']
- # header = nrrd.read_header(com_file)
- # com_position = header['space origin']
- #
- # # with open(file_com) as fobj:
- # # for line in fobj:
- # # line = line.replace('[',']')
- # # line_data = re.split("]|,", line)
- # #
- # # NA = np.asarray(line_data[1:4])
- # # trans_mat = NA.astype(float)
- #
- # # get translation from align from rotation matrix
- # rot_ssm = np.loadtxt(file_align)
- #
- # # rot_mat_ssm = np.transpose(rot_ssm)
- # # rot_mat_ssm = np.vstack((rot_mat_ssm, [0,0,0,1]))
- #
- # # translate points cloud SSM instance to align with original mesh
- # translate = pad_position-com_position+rot_ssm[3,:]
- #
- # # r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output\femur_bone\new_bone\Segmentation_femur_R_short_ssm_reconstruct3.stl' #
- # pre, ext = os.path.splitext(xyz_file)
- # copyfile(particle_file, pre + '.paticles')
- # if not os.path.isfile(pre + '.xyz'):
- # os.rename(pre + '.paticles', pre + '.xyz')
- # mesh3 = xyz_file # =local.particle file
- # ms6 = pymeshlab.MeshSet()
- # ms6.load_new_mesh(mesh3)
- # ms6.apply_filter('transform_translate_center_set_origin', traslmethod =0, axisx=translate[0], axisy=translate[1], axisz=-450)
- # ms6.apply_filter('transform_translate_center_set_origin', traslmethod =0, axisx=0, axisy=0, axisz=-450)
- # ms6.apply_filter('transform_translate_center_set_origin', traslmethod=0, axisx=0, axisy=0, axisz=-450)
- # ms6.apply_filter('transform_translate_center_set_origin', traslmethod =0, axisx=0, axisy=0, axisz=translate[2]+1350)
- # ms6.save_current_mesh(path + '\SSM_' + segment + '_transform.xyz')
+ if short == 1:
+ file_com = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone_short\new_bone\groomed\com_aligned\Segmentation_' + segment + '_' + side + '_short_' + str(
+ subject) + reflect + '.isores.pad.com.txt'
+ file_align = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone_short\new_bone\groomed\aligned\Segmentation_' + segment + '_' + side + '_short_' + str(
+ subject) + reflect + '.isores.pad.com.center.aligned.txt'
+ pad_file = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone_short\new_bone\groomed\padded\segementations\Segmentation_' + segment + '_' + side + '_short_' + str(
+ subject) + reflect + '.isores.pad.nrrd'
+ com_file = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone_short\new_bone\groomed\com_aligned\Segmentation_' + segment + '_' + side + '_short_' + str(
+ subject) + reflect + '.isores.pad.com.nrrd'
+ particle_file = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone_short\new_bone\shape_models\4096\Segmentation_' + segment + '_' + side + '_short_' + str(
+ subject) + reflect + '.isores.pad.com.center.aligned.clipped.cropped.tpSmoothDT_local.particles'
+ xyz_file = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone_short\new_bone\shape_models\Segmentation_' + segment + '_' + side + '_short_' + str(
+ subject) + reflect + '.isores.pad.com.center.aligned.clipped.cropped.tpSmoothDT_local.xyz'
+ # align_file = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output\femur_bone\new_bone\groomed\aligned\Segmentation_femur_' + side + '_short_' + str(subject) + reflect + '.isores.pad.com.center.aligned.nrrd'
+ path_bones = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone_short\new_bone\input'
+ else:
+ file_com = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone\groomed\com_aligned\Segmentation_' + segment + '_' + side + '_short_' + str(subject) + reflect + '.isores.pad.com.txt'
+ file_align = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone\groomed\aligned\Segmentation_' + segment + '_' + side + '_short_' + str(subject) + reflect + '.isores.pad.com.center.aligned.txt'
+ pad_file = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone\groomed\padded\segementations\Segmentation_' + segment + '_' + side + '_short_' + str(subject) + reflect + '.isores.pad.nrrd'
+ com_file = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone\groomed\com_aligned\Segmentation_' + segment + '_' + side + '_short_' + str(subject) + reflect + '.isores.pad.com.nrrd'
+ particle_file = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone\shape_models\4096\Segmentation_' + segment + '_' + side + '_short_' + str(subject) + reflect + '.isores.pad.com.center.aligned.clipped.cropped.tpSmoothDT_local.particles'
+ xyz_file = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone\shape_models\Segmentation_' + segment + '_' + side + '_short_' + str(subject) + reflect + '.isores.pad.com.center.aligned.clipped.cropped.tpSmoothDT_local.xyz'
+ # align_file = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output\femur_bone\new_bone\groomed\aligned\Segmentation_femur_' + side + '_short_' + str(subject) + reflect + '.isores.pad.com.center.aligned.nrrd'
+ path_bones = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone\input'
+
+ # get change in position from nrrd header files
+ header = nrrd.read_header(pad_file)
+ pad_position = header['space origin']
+ header = nrrd.read_header(com_file)
+ com_position = header['space origin']
+
+ # with open(file_com) as fobj:
+ # for line in fobj:
+ # line = line.replace('[',']')
+ # line_data = re.split("]|,", line)
#
- # # run ICP to get final position SSM point cloud on original mesh
+ # NA = np.asarray(line_data[1:4])
+ # trans_mat = NA.astype(float)
+
+ # get translation from align from rotation matrix
+ rot_ssm = np.loadtxt(file_align)
+
+ # rot_mat_ssm = np.transpose(rot_ssm)
+ # rot_mat_ssm = np.vstack((rot_mat_ssm, [0,0,0,1]))
+
+ # translate points cloud SSM instance to align with original mesh
+ translate = pad_position-com_position+rot_ssm[3,:]
+
+ # r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output\femur_bone\new_bone\Segmentation_femur_R_short_ssm_reconstruct3.stl' #
+ pre, ext = os.path.splitext(xyz_file)
+ copyfile(particle_file, pre + '.paticles')
+ if not os.path.isfile(pre + '.xyz'):
+ os.rename(pre + '.paticles', pre + '.xyz')
+ mesh3 = xyz_file # =local.particle file
+ ms6 = pymeshlab.MeshSet()
+ ms6.load_new_mesh(mesh3)
+ ms6.apply_filter('transform_translate_center_set_origin', traslmethod =0, axisx=translate[0], axisy=translate[1], axisz=-430)
+ ms6.apply_filter('transform_translate_center_set_origin', traslmethod =0, axisx=0, axisy=0, axisz=-430)
+ ms6.apply_filter('transform_translate_center_set_origin', traslmethod=0, axisx=0, axisy=0, axisz=-430)
+ ms6.apply_filter('transform_translate_center_set_origin', traslmethod =0, axisx=0, axisy=0, axisz=translate[2]+1290)
+ ms6.save_current_mesh(path + '\SSM_' + segment + '_short_transform.xyz')
+
+ # run ICP to get final position SSM point cloud on original mesh
# mesh = trimesh.load_mesh(path + '\Segmentation_' + segment + '_resample.stl')
- # # mesh = trimesh.load_mesh(path_bones + '\Segmentation_' + segment + '_' + side + '_short_' + str(subject) + '.STL')
- # points = trimesh.load_mesh(path + '\SSM_' + segment + '_transform.xyz')
- # if reflect == '.reflect':
- # M = trimesh.transformations.scale_and_translate((-1,1,1))
- # points.apply_transform(M)
- # kwargs = {"scale": False}
- # icp = trimesh.registration.icp(points.vertices,mesh,initial=np.identity(4),threshold=1e-5,max_iterations=20,**kwargs)
- # points.apply_transform(icp[0])
+ mesh = trimesh.load_mesh(r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output\tibia_bone_short\new_bone\input\Segmentation_' + segment + '_' + side + '_short_' + str(subject) + '_remesh.stl')
+ # mesh = trimesh.load_mesh(path_bones + '\Segmentation_' + segment + '_' + side + '_short_' + str(subject) + '.STL')
+ points = trimesh.load_mesh(path + '\SSM_' + segment + '_short_transform.xyz')
+ if reflect == '.reflect':
+ M = trimesh.transformations.scale_and_translate((-1,1,1))
+ points.apply_transform(M)
+ kwargs = {"scale": False}
+ icp = trimesh.registration.icp(points.vertices,mesh,initial=np.identity(4),threshold=1e-5,max_iterations=20,**kwargs)
+ points.apply_transform(icp[0])
# icp = trimesh.registration.icp(points.vertices, mesh, initial=np.identity(4), threshold=1e-5, max_iterations=20,**kwargs)
# points.apply_transform(icp[0])
- # np.savetxt(path + '\SSM_' + segment + '_transform_icp.xyz', points.vertices, delimiter=" ")
+ np.savetxt(path + '\SSM_' + segment + '_short_transform_icp.xyz', points.vertices, delimiter=" ")
# ms5 = pymeshlab.MeshSet()
# ms5.load_new_mesh(path + '\SSM_' + segment + '_transform_icp.xyz')
- points = trimesh.load_mesh(path + '\SSM_' + segment + '_transform_icp.xyz')
+ points = trimesh.load_mesh(path + '\SSM_' + segment + '_short_transform_icp.xyz')
Counter = len(glob.glob1(path, 'Segmentation_' + segment + '_area*.stl'))
close_verts = []
close_verts_verts = np.empty([0,3])
@@ -123,10 +148,10 @@ for segment in segments:
# close_verts = np.vstack([close_verts, np.where(abs(distance)<2)])
# close_verts_verts.append(points.vertices[np.where(abs(distance)<2)])
close_verts_verts = np.vstack([close_verts_verts, points.vertices[np.where(abs(distance)<1)]])
- np.savetxt(path + '\SSM_' + segment + '_areas.xyz', np.asarray(close_verts_verts), delimiter=" ")
+ np.savetxt(path + '\SSM_' + segment + '_short_areas.xyz', np.asarray(close_verts_verts), delimiter=" ")
- for lig in range(0, 8):
+ for lig in range(0, 11):
lig_no = ligaments[lig][ind]
if not lig_no == 0:
SSMpoints[lig][ind] = close_verts[lig_no-1][0]
@@ -135,28 +160,48 @@ for segment in segments:
from collections import Counter
occurances = []
all_occ = []
- for ind in range(0,7):
+ orders = []
+ for ind in range(0,11):
count=0
- occur = Counter(np.concatenate(SSMpoints[ind], axis=0))
+ occur = []
+ if ind == 2:
+ occur = Counter(np.concatenate(SSMpoints[ind]+SSMpoints[ind+1], axis=0))
+ elif ind == 4:
+ occur = Counter(np.concatenate(SSMpoints[ind]+ SSMpoints[ind + 1]+ SSMpoints[ind + 2]+ SSMpoints[ind + 3], axis=0))
+ elif ind == 3 or ind == 5 or ind == 6 or ind == 7:
+ continue
+ else:
+ occur = Counter(np.concatenate(SSMpoints[ind], axis=0))
order = occur.most_common()
for j in range(0,10):
if len(SSMpoints[ind][j]):
count=count+1
+ if ind == 4:
+ print(count)
try:
- index = [x for x, y in enumerate(order) if y[1] == count/2]
- most_occ = order[0:index[0]-1]
+ index = [x for x, y in enumerate(order) if y[1] == int(count/2)]
+ most_occ = order[0:index[-1]]
except:
most_occ = order[0:1]
all_occ.append(np.asarray([i[0] for i in order]))
occurances.append(np.asarray([i[0] for i in most_occ]))
+ orders.append(np.asarray([i[1] for i in order]))
- points = trimesh.load_mesh(r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone\shape_models\mean_shape.xyz')
- np.savetxt(r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone\shape_models\meanshape_ligs.xyz',
+ points = trimesh.load_mesh(r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone_short\new_bone\shape_models\mean_shape.xyz')
+ np.savetxt(r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone_short\new_bone\shape_models\meanshape_ligs.xyz',
points.vertices[np.hstack(occurances).astype(int)], delimiter=" ")
+
+ pred_lig_points_color = np.c_[points.vertices[np.hstack(all_occ).astype(int)], np.hstack(orders).astype(int)]
+ np.savetxt(
+ r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone_short\new_bone\shape_models\meanshape_ligs_color.xyz',
+ pred_lig_points_color, delimiter=" ")
+
mask = np.ones(len(points.vertices), dtype=bool)
mask[np.hstack(occurances).astype(int)] = False
- np.savetxt(r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone\shape_models\meanshape_bone_no_lig.xyz',
+ np.savetxt(r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone_short\new_bone\shape_models\meanshape_bone_no_lig.xyz',
points.vertices[np.where(mask)], delimiter=" ")
+ np.savez(r'C:\Users\mariskawesseli\Documents\LigamentStudy\ImageData\occurances_order_short.npy', PCL=orders[0],MCLp=orders[1],
+ MCLd=orders[2],post_obl=orders[3],ACL=orders[4],LCL=orders[5],pop=orders[6])
# print('Surface area femur ligament' + str(lig_no) + ': ' + str(surface) + ' mm2')
@@ -176,9 +221,12 @@ for segment in segments:
side = 'L'
reflect = '.reflect'
- points = trimesh.load_mesh(path + '\SSM_' + segment + '_transform_icp.xyz')
+ points = trimesh.load_mesh(path + '\SSM_' + segment + '_short_transform_icp.xyz')
pred_lig_points = points.vertices[np.hstack(occurances).astype(int)]
- np.savetxt(path + '\SSM_' + segment + '_pred_points.xyz', np.asarray(pred_lig_points), delimiter=" ")
+ np.savetxt(path + '\SSM_' + segment + '_short_pred_points.xyz', np.asarray(pred_lig_points), delimiter=" ")
+ pred_lig_points_color = np.c_[points.vertices[np.hstack(all_occ).astype(int)],np.hstack(orders).astype(int)]
+ np.savetxt(path + '\SSM_' + segment + '_short_pred_points_color.xyz', np.asarray(pred_lig_points_color), delimiter=" ")
+
- np.savez(r'C:\Users\mariskawesseli\Documents\LigamentStudy\ImageData\occurances_' + segment + '.npy',PCL=occurances[0],MCLp=occurances[1],
+ np.savez(r'C:\Users\mariskawesseli\Documents\LigamentStudy\ImageData\occurances_' + segment + '_short',PCL=occurances[0],MCLp=occurances[1],
MCLd=occurances[2],post_obl=occurances[3],ACL=occurances[4],LCL=occurances[5],pop=occurances[6])
diff --git a/LigamentStudy/fitSSM_mri.py b/LigamentStudy/fitSSM_mri.py
new file mode 100644
index 0000000..3a6eb4d
--- /dev/null
+++ b/LigamentStudy/fitSSM_mri.py
@@ -0,0 +1,219 @@
+import pymeshlab
+import numpy as np
+import trimesh
+import nrrd
+import re
+import os
+import pandas as pd
+from tabulate import tabulate
+from shutil import copyfile
+import glob
+
+subjects = ['S0'] # [9,13,19,23,26,29,32,35,37,41]
+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],
+ [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]]
+
+for segment in segments:
+ if segment == 'femur':
+ ligaments = ligaments_fem
+ else:
+ ligaments = ligaments_tib
+
+ SSMpoints = [[] for i in range(11)]
+ for ind in range(0,11):
+ SSMpoints[ind] = [[] for i in range(10)]
+
+ for ind, subject in enumerate(subjects):
+ path = os.path.join(r"C:\Users\mariskawesseli\Documents\LigamentStudy\MRI\S0_prelim")
+ if subject in [9,13,26,29,32]:
+ side = 'R'
+ reflect = ''
+ else:
+ side = 'L'
+ reflect = '.reflect'
+
+ """SSM part"""
+ # files from SSM workflow shapeworks
+
+ file_com = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone_short\new_bone_mri\groomed\com_aligned\S0_' + side + '_' + segment + reflect + '.isores.pad.com.txt'
+ file_align = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone_short\new_bone_mri\groomed\aligned\S0_' + side + '_' + segment + reflect + '.isores.pad.com.center.aligned.txt'
+ pad_file = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone_short\new_bone_mri\groomed\padded\segementations\S0_' + side + '_' + segment + reflect + '.isores.pad.nrrd'
+ com_file = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone_short\new_bone_mri\groomed\com_aligned\S0_' + side + '_' + segment + reflect + '.isores.pad.com.nrrd'
+ if segment == 'femur':
+ particle_file = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone_mri\shape_models\4096\S0_' + side + '_' + segment + reflect + '.isores.pad.com.center.aligned.clipped.cropped.tpSmoothDT_local.particles'
+ else:
+ particle_file = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone_short\new_bone_mri\shape_models\4096\S0_' + side + '_' + segment + reflect + '.isores.pad.com.center.aligned.clipped.cropped.tpSmoothDT_local.particles'
+ xyz_file = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone_short\new_bone_mri\shape_models\S0_' + side + '_' + segment + reflect + '.isores.pad.com.center.aligned.clipped.cropped.tpSmoothDT_local.xyz'
+ # align_file = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output\femur_bone\new_bone\groomed\aligned\Segmentation_femur_' + side + '_short_' + str(subject) + reflect + '.isores.pad.com.center.aligned.nrrd'
+ path_bones = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone_short\new_bone_mri\input'
+
+ # get change in position from nrrd header files
+ header = nrrd.read_header(pad_file)
+ pad_position = header['space origin']
+ header = nrrd.read_header(com_file)
+ com_position = header['space origin']
+
+ # with open(file_com) as fobj:
+ # for line in fobj:
+ # line = line.replace('[',']')
+ # line_data = re.split("]|,", line)
+ #
+ # NA = np.asarray(line_data[1:4])
+ # trans_mat = NA.astype(float)
+
+ # get translation from align from rotation matrix
+ rot_ssm = np.loadtxt(file_align)
+
+ # rot_mat_ssm = np.transpose(rot_ssm)
+ # rot_mat_ssm = np.vstack((rot_mat_ssm, [0,0,0,1]))
+
+ # translate points cloud SSM instance to align with original mesh
+ translate = pad_position-com_position+rot_ssm[3,:]
+
+ # r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output\femur_bone\new_bone\Segmentation_femur_R_short_ssm_reconstruct3.stl' #
+ pre, ext = os.path.splitext(xyz_file)
+ copyfile(particle_file, pre + '.paticles')
+
+ if not os.path.isfile(pre + '.xyz'):
+ os.rename(pre + '.paticles', pre + '.xyz')
+ mesh3 = xyz_file # =local.particle file
+ ms6 = pymeshlab.MeshSet()
+ ms6.load_new_mesh(mesh3)
+ ms6.apply_filter('transform_translate_center_set_origin', traslmethod=0, axisx=translate[0], axisy=translate[1], axisz=translate[2])
+ ms6.save_current_mesh(path + '\SSM_' + segment + '_transform.xyz')
+
+ # run ICP to get final position SSM point cloud on original mesh
+ mesh = trimesh.load_mesh(path + '/bone_tibia_2_bone_rot.stl') # bone_femur2_2_bone_rot.stl
+ # mesh = trimesh.load_mesh(path_bones + '\Segmentation_' + segment + '_' + side + '_short_' + str(subject) + '.STL')
+ points = trimesh.load_mesh(path + '\SSM_' + segment + '_transform.xyz')
+ if reflect == '.reflect':
+ M = trimesh.transformations.scale_and_translate((-1,1,1))
+ points.apply_transform(M)
+ kwargs = {"scale": False}
+ icp = trimesh.registration.icp(points.vertices,mesh,initial=np.identity(4),threshold=1e-5,max_iterations=20,**kwargs)
+ points.apply_transform(icp[0])
+ icp = trimesh.registration.icp(points.vertices, mesh, initial=np.identity(4), threshold=1e-5, max_iterations=20,**kwargs)
+ points.apply_transform(icp[0])
+ np.savetxt(path + '\SSM_' + segment + '_transform_icp.xyz', points.vertices, delimiter=" ")
+
+ # ms5 = pymeshlab.MeshSet()
+ # ms5.load_new_mesh(path + '\SSM_' + segment + '_transform_icp.xyz')
+ # points = trimesh.load_mesh(path + '\SSM_' + segment + '_transform_icp.xyz')
+ # Counter = len(glob.glob1(path, 'Segmentation_' + segment + '_area*.stl'))
+ # close_verts = []
+ # close_verts_verts = np.empty([0,3])
+ # for count in range(1, Counter + 1):
+ # mesh = trimesh.load_mesh(os.path.join(path,'Segmentation_' + segment + '_area' + str(count) + '.stl'))
+ # # [closest, distance, id] = trimesh.proximity.closest_point(mesh, points.vertices)
+ # distance = trimesh.proximity.signed_distance(mesh, points.vertices)
+ # close_verts.append(np.where(abs(distance)<1))
+ # # close_verts = np.vstack([close_verts, np.where(abs(distance)<2)])
+ # # close_verts_verts.append(points.vertices[np.where(abs(distance)<2)])
+ # close_verts_verts = np.vstack([close_verts_verts, points.vertices[np.where(abs(distance)<1)]])
+ # np.savetxt(path + '\SSM_' + segment + '_areas.xyz', np.asarray(close_verts_verts), delimiter=" ")
+
+
+ # for lig in range(0, 11):
+ # lig_no = ligaments[lig][ind]
+ # if not lig_no == 0:
+ # SSMpoints[lig][ind] = close_verts[lig_no-1][0]
+ #
+ # # dupes = [x for n, x in enumerate(np.concatenate(SSMpoints[0])) if x in np.concatenate(SSMpoints[0])[:n]]
+ # from collections import Counter
+ # occurances = []
+ # all_occ = []
+ # orders = []
+ # for ind in range(0,11):
+ # count=0
+ # occur = []
+ # if ind == 2:
+ # occur = Counter(np.concatenate(SSMpoints[ind]+SSMpoints[ind+1], axis=0))
+ # elif ind == 4:
+ # occur = Counter(np.concatenate(SSMpoints[ind]+ SSMpoints[ind + 1]+ SSMpoints[ind + 2]+ SSMpoints[ind + 3], axis=0))
+ # elif ind == 3 or ind == 5 or ind == 6 or ind == 7:
+ # continue
+ # else:
+ # occur = Counter(np.concatenate(SSMpoints[ind], axis=0))
+ # order = occur.most_common()
+ # for j in range(0,10):
+ # if len(SSMpoints[ind][j]):
+ # count=count+1
+ # try:
+ # index = [x for x, y in enumerate(order) if y[1] == int(count/2)]
+ # most_occ = order[0:index[0]-1]
+ # except:
+ # most_occ = order[0:1]
+ # all_occ.append(np.asarray([i[0] for i in order]))
+ # occurances.append(np.asarray([i[0] for i in most_occ]))
+ # orders.append(np.asarray([i[1] for i in order]))
+ #
+ # points = trimesh.load_mesh(r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone\shape_models\mean_shape.xyz')
+ # np.savetxt(r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone\shape_models\meanshape_ligs.xyz',
+ # points.vertices[np.hstack(occurances).astype(int)], delimiter=" ")
+ #
+ # pred_lig_points_color = np.c_[points.vertices[np.hstack(all_occ).astype(int)], np.hstack(orders).astype(int)]
+ # np.savetxt(
+ # r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone\shape_models\meanshape_ligs_color.xyz',
+ # pred_lig_points_color, delimiter=" ")
+ #
+ # mask = np.ones(len(points.vertices), dtype=bool)
+ # mask[np.hstack(occurances).astype(int)] = False
+ # np.savetxt(r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone\shape_models\meanshape_bone_no_lig.xyz',
+ # points.vertices[np.where(mask)], delimiter=" ")
+ # np.savez(r'C:\Users\mariskawesseli\Documents\LigamentStudy\ImageData\occurances_order.npy', PCL=orders[0],MCLp=orders[1],
+ # MCLd=orders[2],post_obl=orders[3],ACL=orders[4],LCL=orders[5],pop=orders[6])
+ #
+ #
+ # # print('Surface area femur ligament' + str(lig_no) + ': ' + str(surface) + ' mm2')
+ # # ms5.load_new_mesh(os.path.join(path,'Segmentation_femur_area' + str(count) + '.stl'))
+ # # no_meshes = ms5.number_meshes()
+ # # ms5.apply_filter('distance_from_reference_mesh', measuremesh=0, refmesh=no_meshes-1, signeddist=False)
+ # # ms5.set_current_mesh(0)
+ # # ms5.conditional_vertex_selection(condselect="q<1")
+ # # m = ms5.current_mesh()
+
+ occ = np.load(r'C:\Users\mariskawesseli\Documents\LigamentStudy\ImageData\occurances_' + segment + '_short.npz')
+ occurances = [occ['PCL'], occ['MCLp'], occ['MCLd'], occ['post_obl'], occ['ACL'], occ['LCL'], occ['pop']]
+
+ subjects = ['S0'] # [9,13,19,23,26,29,32,35,37,41]
+ for ind, subject in enumerate(subjects):
+ path = os.path.join(r"C:\Users\mariskawesseli\Documents\LigamentStudy\MRI\S0_prelim")
+ if subject in [9,13,26,29,32]:
+ side = 'R'
+ reflect = ''
+ else:
+ side = 'L'
+ reflect = '.reflect'
+
+ points = trimesh.load_mesh(path + '\SSM_' + segment + '_transform_icp.xyz')
+ pred_lig_points = points.vertices[np.hstack(occurances).astype(int)]
+ np.savetxt(path + '\SSM_' + segment + '_short_pred_points.xyz', np.asarray(pred_lig_points), delimiter=" ")
+ # pred_lig_points_color = np.c_[points.vertices[np.hstack(all_occ).astype(int)],np.hstack(orders).astype(int)]
+ # np.savetxt(path + '\SSM_' + segment + '_pred_points_color.xyz', np.asarray(pred_lig_points_color), delimiter=" ")
+
+
+ # np.savez(r'C:\Users\mariskawesseli\Documents\LigamentStudy\ImageData\occurances_' + segment + '.npy',PCL=occurances[0],MCLp=occurances[1],
+ # MCLd=occurances[2],post_obl=occurances[3],ACL=occurances[4],LCL=occurances[5],pop=occurances[6])
+
diff --git a/LigamentStudy/scaleOsim.py b/LigamentStudy/scaleOsim.py
new file mode 100644
index 0000000..9bbeaba
--- /dev/null
+++ b/LigamentStudy/scaleOsim.py
@@ -0,0 +1,170 @@
+import trimesh
+import numpy as np
+from scipy import interpolate
+import matplotlib.pyplot as plt
+
+def interpolate_lig_points(lig_points, no_points, plot = 1):
+ tck, u = interpolate.splprep([lig_points[:, 0], lig_points[:, 1], lig_points[:, 2]],s=10,k=1)
+ x_knots, y_knots, z_knots = interpolate.splev(tck[0], tck)
+ u_fine = np.linspace(0, 1, no_points)
+ x_fine, y_fine, z_fine = interpolate.splev(u_fine, tck, der=0)
+ lig_points_osim = np.transpose(np.asarray([x_fine, y_fine, z_fine]))
+
+ if plot == 1:
+ fig2 = plt.figure()
+ ax3d = fig2.add_subplot(111, projection='3d')
+ ax3d.plot(lig_points[:, 0], lig_points[:, 1], lig_points[:, 2], 'r*')
+ ax3d.plot(x_knots, y_knots, z_knots, 'bo')
+ ax3d.plot(x_fine, y_fine, z_fine, 'go')
+ fig2.show()
+ plt.show()
+
+ return lig_points_osim
+
+
+"""femur"""
+# # run ICP to get final position SSM point cloud on original mesh
+# mesh1 = trimesh.load_mesh('C:\opensim-jam\jam-resources\jam-resources-main\models\knee_healthy\smith2019\Geometry\smith2019-R-femur-bone_remesh.stl')
+# mesh2 = trimesh.load_mesh(r'C:\Users\mariskawesseli\Documents\LigamentStudy\MRI\S0_prelim\bone_femur2_2_bone_rot_remesh2.STL')
+# # points = trimesh.load_mesh(path + '\SSM_' + segment + '_transform.xyz')
+#
+# M = trimesh.transformations.scale_and_translate((1,1,-1))
+# mesh1.apply_transform(M)
+# T = trimesh.transformations.translation_matrix([64.724205, -26.297621, -95.929390])
+# origin, xaxis, yaxis, zaxis = [0,0,0], [1, 0, 0], [0, 1, 0], [0, 0, 1]
+# Rx = trimesh.transformations.rotation_matrix(-90/(180/np.pi), xaxis)
+# Ry = trimesh.transformations.rotation_matrix(90/(180/np.pi), yaxis)
+# R = trimesh.transformations.concatenate_matrices(T, Ry, Rx)
+# mesh2.apply_transform(R)
+# mesh2.export(r'C:\Users\mariskawesseli\Documents\LigamentStudy\MRI\S0_prelim\bone_femur2_2_bone_rot_test.STL')
+#
+# s, fi = trimesh.sample.sample_surface_even(mesh2, 32015, radius=None)
+# mesh3 = trimesh.icp[0]imesh(vertices=s, faces=mesh2.faces[fi])
+# mesh3.export(r'C:\Users\mariskawesseli\Documents\LigamentStudy\MRI\S0_prelim\bone_femur2_2_bone_rot_test.STL')
+# sort_index = np.argsort(fi) # sort and group the face indices
+# points = s[sort_index, :]
+# faceIndices = fi[sort_index]
+# uniqueFaceIndices = np.unique(faceIndices)
+# allMeshPatches = trimesh.icp[0]imesh()
+# pointGroups = [points[faceIndices == i] for i in uniqueFaceIndices]
+# for faceIndex, pointsOnFace in zip(uniqueFaceIndices, pointGroups):
+# meshpatch = trimesh.icp[0]imesh(mesh2.vertices[mesh2.faces[faceIndex, :]].reshape(3, 3),
+# np.array([0, 1, 2]).reshape(1, 3))
+# allMeshPatches += meshpatch
+# allMeshPatches.export(r'C:\Users\mariskawesseli\Documents\LigamentStudy\MRI\S0_prelim\bone_femur2_2_bone_rot_remesh_test.STL')
+#
+# kwargs = {"scale": False}
+# icp1 = trimesh.registration.icp(mesh2.vertices,mesh1,initial=np.identity(4),threshold=1e-5,max_iterations=20,**kwargs)
+#
+# kwargs = {"scale": True}
+# mesh2.export(r'C:\Users\mariskawesseli\Documents\LigamentStudy\MRI\S0_prelim\bone_femur2_2_bone_rot_icp.STL')
+# icp = trimesh.registration.icp(mesh2.vertices, mesh1, initial=icp1[0], threshold=1e-5, max_iterations=20,**kwargs)
+# mesh2.apply_transform(icp[0])
+# mesh2.export(r'C:\Users\mariskawesseli\Documents\LigamentStudy\MRI\S0_prelim\bone_femur2_2_bone_rot_icp.STL')
+# scale, shear, angles, trans, persp = trimesh.transformations.decompose_matrix(icp[0])
+# # mesh1.export(r'C:\Users\mariskawesseli\Documents\LigamentStudy\MRI\S0_prelim\test.STL')
+#
+# ligs = trimesh.load_mesh(r'C:\Users\mariskawesseli\Documents\LigamentStudy\MRI\S0_prelim\femur\SSM_femur_pred_points.xyz')
+# ligs.apply_transform(R)
+# ligs.apply_transform(icp[0])
+# np.savetxt(r'C:\Users\mariskawesseli\Documents\LigamentStudy\MRI\S0_prelim\femur\SSM_femur_pred_points_osim.xyz', ligs.vertices, delimiter=" ")
+
+
+ligs = trimesh.load_mesh(r'C:\Users\mariskawesseli\Documents\LigamentStudy\MRI\S0_prelim\femur\SSM_femur_pred_points_osim.xyz')
+PCL = ligs.vertices[0:61]
+al = [0,1,2,3,4,5,7,10,11,13,29,30,31,34,35,39,40,42,43,47,49,50,51,55,56,57,58,59,24,27,28,48]
+pm = [item for item in list(range(61)) if item not in al]
+PCLal_osim = interpolate_lig_points(PCL[al,:],5)
+PCLpm_osim = interpolate_lig_points(PCL[pm,:],5)
+MCLs = ligs.vertices[61:71]
+MCLs_osim = interpolate_lig_points(MCLs,6)
+MCLd = ligs.vertices[71:73]
+MCLd_osim = interpolate_lig_points(MCLd,5)
+post_obl = ligs.vertices[73:81]
+post_obl_osim = interpolate_lig_points(post_obl,5)
+ACL = ligs.vertices[81:100]
+al = [0,1,2,4,7,9,12,15,18]
+pm = [item for item in list(range(19)) if item not in al]
+ACLal_osim = interpolate_lig_points(ACL[al,:],6)
+ACLpm_osim = interpolate_lig_points(ACL[pm,:],6)
+LCL = ligs.vertices[100:105]
+LCL_osim = interpolate_lig_points(LCL,4)
+
+osim_points = np.concatenate([np.asarray(PCLal_osim),np.asarray(PCLpm_osim), np.asarray(MCLs_osim),
+ np.asarray(MCLd_osim), np.asarray(post_obl_osim), np.asarray(ACLal_osim),
+ np.asarray(ACLpm_osim), np.asarray(LCL_osim)])
+np.savetxt(r'C:\Users\mariskawesseli\Documents\LigamentStudy\MRI\S0_prelim\femur\SSM_femur_pred_points_osim_interp.xyz', osim_points, delimiter=" ")
+
+# fig2 = plt.figure()
+# ax3d = fig2.add_subplot(111, projection='3d')
+# ax3d.plot(ACL[al, 0], ACL[al, 1], ACL[al, 2], 'r*')
+# ax3d.plot(ACL[pm, 0], ACL[pm, 1], ACL[pm, 2], 'bo')
+# fig2.show()
+# plt.show()
+
+"""tibia"""
+# # run ICP to get final position SSM point cloud on original mesh
+# mesh1 = trimesh.load_mesh('C:\opensim-jam\jam-resources\jam-resources-main\models\knee_healthy\smith2019\Geometry\smith2019-R-tibia-bone_remesh.stl')
+# mesh2 = trimesh.load_mesh(r'C:\Users\mariskawesseli\Documents\LigamentStudy\MRI\S0_prelim\bone_tibia_2_bone_rot_remesh.STL')
+# # points = trimesh.load_mesh(path + '\SSM_' + segment + '_transform.xyz')
+#
+# M = trimesh.transformations.scale_and_translate((1,1,-1))
+# mesh1.apply_transform(M)
+# T = trimesh.transformations.translation_matrix([101.562462, -72.768566, -17.893391])
+# origin, xaxis, yaxis, zaxis = [0,0,0], [1, 0, 0], [0, 1, 0], [0, 0, 1]
+# Rx = trimesh.transformations.rotation_matrix(-90/(180/np.pi), xaxis)
+# Ry = trimesh.transformations.rotation_matrix(90/(180/np.pi), yaxis)
+# R = trimesh.transformations.concatenate_matrices(Ry, Rx)
+# mesh2.apply_transform(T)
+# mesh2.apply_transform(R)
+# mesh2.export(r'C:\Users\mariskawesseli\Documents\LigamentStudy\MRI\S0_prelim\bone_tibia_2_bone_rot_remesh_test.STL')
+#
+# kwargs = {"scale": False}
+# icp1 = trimesh.registration.icp(mesh2.vertices,mesh1,initial=np.identity(4),threshold=1e-5,max_iterations=20,**kwargs)
+#
+# kwargs = {"scale": True}
+# # mesh2.export(r'C:\Users\mariskawesseli\Documents\LigamentStudy\MRI\S0_prelim\bone_tibia_2_bone_rot_icp.STL')
+# icp = trimesh.registration.icp(mesh2.vertices, mesh1, initial=icp1[0], threshold=1e-5, max_iterations=20,**kwargs)
+# mesh2.apply_transform(icp[0])
+# mesh2.export(r'C:\Users\mariskawesseli\Documents\LigamentStudy\MRI\S0_prelim\bone_tibia_2_bone_rot_icp.STL')
+# scale, shear, angles, trans, persp = trimesh.transformations.decompose_matrix(icp[0])
+# mesh1.export(r'C:\Users\mariskawesseli\Documents\LigamentStudy\MRI\S0_prelim\test_tib.STL')
+#
+# ligs = trimesh.load_mesh(r'C:\Users\mariskawesseli\Documents\LigamentStudy\MRI\S0_prelim\SSM_tibia_short_pred_points.xyz')
+# ligs.apply_transform(T)
+# ligs.apply_transform(R)
+# ligs.apply_transform(icp[0])
+# np.savetxt(r'C:\Users\mariskawesseli\Documents\LigamentStudy\MRI\S0_prelim\tibia\SSM_tibia_short_pred_points_osim.xyz', ligs.vertices, delimiter=" ")
+
+
+# ligs = trimesh.load_mesh(r'C:\Users\mariskawesseli\Documents\LigamentStudy\MRI\S0_prelim\tibia\SSM_tibia_short_pred_points_osim.xyz')
+# PCL = ligs.vertices[0:71]
+# al = [0,1,2,4,6,7,10,12,13,19,21,22,23,24,25,26,27,31,32,34,35,40,45,46,47,50,52,56,59,60,61,62,63,65,70]
+# pm = [item for item in list(range(71)) if item not in al]
+# PCLal_osim = interpolate_lig_points(PCL[al,:],5)
+# PCLpm_osim = interpolate_lig_points(PCL[pm,:],5)
+# # MCLs = ligs.vertices[61:72]
+# # MCLs_osim = interpolate_lig_points(MCLs,6)
+# MCLd = ligs.vertices[71:82]
+# MCLd_osim = interpolate_lig_points(MCLd,5)
+# post_obl = ligs.vertices[82:86]
+# post_obl_osim = interpolate_lig_points(post_obl,5)
+# ACL = ligs.vertices[86:150]
+# al = [0,1,2,4,6,7,10,12,13,19,21,22,23,24,25,26,27,31,32,34,35,40,45,46,47,50,52,56,59,60,61,62,63]
+# pm = [item for item in list(range(64)) if item not in al]
+# ACLal_osim = interpolate_lig_points(ACL[al,:],6)
+# ACLpm_osim = interpolate_lig_points(ACL[pm,:],6)
+# # LCL = ligs.vertices[101:106]
+# # LCL_osim = interpolate_lig_points(LCL,4)
+#
+# osim_points = np.concatenate([np.asarray(PCLal_osim),np.asarray(PCLpm_osim),
+# np.asarray(MCLd_osim), np.asarray(post_obl_osim), np.asarray(ACLal_osim),
+# np.asarray(ACLpm_osim)])
+# np.savetxt(r'C:\Users\mariskawesseli\Documents\LigamentStudy\MRI\S0_prelim\tibia\SSM_tibia_pred_points_osim_interp.xyz', osim_points, delimiter=" ")
+
+# fig2 = plt.figure()
+# ax3d = fig2.add_subplot(111, projection='3d')
+# ax3d.plot(ligs.vertices[:, 0], ligs.vertices[:, 1], ligs.vertices[:, 2], 'r*')
+# ax3d.plot(LCL[:, 0], LCL[:, 1], LCL[:, 2], 'bo')
+# fig2.show()
+# plt.show()
\ No newline at end of file
--
GitLab