From e27937e1cdc2a00dfcf628f02d4864e37e196bf3 Mon Sep 17 00:00:00 2001
From: Mariska Wesseling <m.g.h.wesseling@tudelft.nl>
Date: Thu, 27 Jan 2022 14:03:02 +0100
Subject: [PATCH] cleaned up fitSSM_mri.py
---
LigamentStudy/Analyses.py | 67 +++++-----
LigamentStudy/ReadMe.md | 27 ++++
LigamentStudy/fitSSM_mri.py | 253 ++++++------------------------------
LigamentStudy/scaleOsim.py | 133 +++++++++----------
4 files changed, 157 insertions(+), 323 deletions(-)
create mode 100644 LigamentStudy/ReadMe.md
diff --git a/LigamentStudy/Analyses.py b/LigamentStudy/Analyses.py
index c0f60d2..a6dc900 100644
--- a/LigamentStudy/Analyses.py
+++ b/LigamentStudy/Analyses.py
@@ -31,41 +31,42 @@ from openpyxl import load_workbook
subjects = [9,13,19,23,26,29,32,35,37,41] #9,13,19,23,26,29,32,35,41
segments = ['tibia','femur'] #'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],
- [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]]
+# the mesh number related to the ligament for each specimen
+ligaments_fem = [[1,1,1,1,1,1,1,1,1,1], # PCL
+ [6,5,6,6,6,6,4,4,5,5], # MCLp
+ [3,2,5,3,3,2,2,0,3,3], # MCLd
+ [0,8,0,0,0,0,0,0,0,0], # MCLd2
+ [7,3,7,7,7,5,7,6,7,0], # POL
+ [0,0,8,0,0,0,0,0,0,0], # POL2
+ [0,0,0,0,0,0,0,0,0,0], # POL3
+ [0,0,0,0,0,0,0,0,0,0], # POL4
+ [4,6,3,5,4,0,0,3,4,4], # ACL
+ [5,7,4,4,5,7,6,5,6,6], # LCL
+ [2,4,2,2,2,3,3,2,2,2]] # POP
-ligaments_tib = [[5, 7, 6, 5, 3, 4, 4, 5, 5, 4],
- [1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
- [3, 3, 8, 3, 5, 3, 5, 0, 3, 3],
- [0, 4, 0, 0, 0, 0, 0, 0, 0, 0], # MCLd2
- [4, 5, 3, 4, 4, 5, 3, 2, 4, 0],
- [0, 6, 4, 0, 0, 0, 0, 0, 0, 0], # POL2
- [0, 0, 5, 0, 0, 0, 0, 0, 0, 0], # POL3
- [0, 0, 7, 0, 0, 0, 0, 0, 0, 0], # POL4
- [6, 8, 9, 6, 6, 6, 6, 6, 6, 5],
- [2, 2, 2, 2, 2, 2, 2, 3, 2, 2],
- [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]
+ligaments_tib = [[5,7,6,5,3,4,4,5,5,4], # PCL
+ [1,1,1,1,1,1,1,1,1,1], # MCLp
+ [3,3,8,3,5,3,5,0,3,3], # MCLd
+ [0,4,0,0,0,0,0,0,0,0], # MCLd2
+ [4,5,3,4,4,5,3,2,4,0], # POL
+ [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], # ACL
+ [2,2,2,2,2,2,2,3,2,2], # LCL
+ [0,0,0,0,0,0,0,0,0,0]] # POP
-ligaments_fib = [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0], # PCL
- [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], # MCLp
- [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], # MCLd
- [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], # MCLd2
- [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], # POL
- [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], # POL2
- [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], # POL3
- [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], # POL4
- [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], # ACL
- [2, 2, 2, 2, 2, 2, 2, 3, 2, 2], # LCL
- [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]] # POP
+ligaments_fib = [[0,0,0,0,0,0,0,0,0,0], # PCL
+ [0,0,0,0,0,0,0,0,0,0], # MCLp
+ [0,0,0,0,0,0,0,0,0,0], # MCLd
+ [0,0,0,0,0,0,0,0,0,0], # MCLd2
+ [0,0,0,0,0,0,0,0,0,0], # POL
+ [0,0,0,0,0,0,0,0,0,0], # POL2
+ [0,0,0,0,0,0,0,0,0,0], # POL3
+ [0,0,0,0,0,0,0,0,0,0], # POL4
+ [0,0,0,0,0,0,0,0,0,0], # ACL
+ [2,2,2,2,2,2,2,3,2,2], # LCL
+ [0,0,0,0,0,0,0,0,0,0]] # POP
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')
diff --git a/LigamentStudy/ReadMe.md b/LigamentStudy/ReadMe.md
new file mode 100644
index 0000000..635bec0
--- /dev/null
+++ b/LigamentStudy/ReadMe.md
@@ -0,0 +1,27 @@
+
+X. Fit SSM to new bone
+Make sure you have a SSM of the specific bone in Shapeworks and know which points represent the ligament attachments.
+
+1. Segment bone from images
+2. Groom pipeline for new bone
+Reflect bone if left (put _L in the name) as SSM is for right bones.
+Make sure the bone is correctly aligned and cropped as the input bones of the original SSM.
+For cropping use the median bone of the original SSM
+If needed adapt the padding
+3. Create the mean shape particle file of the original SSM (getMeanShape.py)
+4. Run the SSM with fixed domains, where the number of fixed domains is the number of original input bones used for the SSM.
+If needed increase the narrow band
+5. Position the particle file (SSM point cloud) at the original bone location (FitSSM_mri.py)
+Make sure the original bone is not too large, to avoid very slow ICP. Aim for a file size below 10MB (about 20,000 faces) (Meshlab - Simplificaion: Quadric Edge Collapse Decimation)
+Required variables:
+subjects => name of the subjects you want to process
+sides => for each subject the side that is being analyzed
+segments => segments you want to analyze
+short_ssm => for each segment, is the shorter SSM needed (0=false, 1=true)
+no_particles => for each segment, number of particles in the SSM
+6. Get the points associated with all ligament locations on the original bone mesh location
+7. For each ligament, determine the SSM points associated to the ligament attachments (scaleOsim.py and adaptLigaments.py)
+Interpolate the points to obtain the number of points needed for the OpenSim model
+Write points to osim file
+8. Scale the ligament parameter based on the length of the ligament (still in Matlab)
+
diff --git a/LigamentStudy/fitSSM_mri.py b/LigamentStudy/fitSSM_mri.py
index db7c61e..e00ebe7 100644
--- a/LigamentStudy/fitSSM_mri.py
+++ b/LigamentStudy/fitSSM_mri.py
@@ -10,242 +10,65 @@ from shutil import copyfile
import glob
subjects = ['S0'] # [9,13,19,23,26,29,32,35,37,41]
-segments = ['fibula'] #'femur',
-ligaments_fem = [[1,1,1,1,1,1,1,1,1,1], # PCL
- [6,5,6,6,6,6,4,4,5,5], # MCLp
- [3,2,5,3,3,2,2,0,3,3], # MCLd
- [0,8,0,0,0,0,0,0,0,0], # MCLd2
- [7,3,7,7,7,5,7,6,7,0], # POL
- [0,0,8,0,0,0,0,0,0,0], # POL2
- [0,0,0,0,0,0,0,0,0,0], # POL3
- [0,0,0,0,0,0,0,0,0,0], # POL4
- [4,6,3,5,4,0,0,3,4,4], # ACL
- [5,7,4,4,5,7,6,5,6,6], # LCL
- [2,4,2,2,2,3,3,2,2,2]] # POP
-
-ligaments_tib = [[5,7,6,5,3,4,4,5,5,4], # PCL
- [1,1,1,1,1,1,1,1,1,1], # MCLp
- [3,3,8,3,5,3,5,0,3,3], # MCLd
- [0,4,0,0,0,0,0,0,0,0], # MCLd2
- [4,5,3,4,4,5,3,2,4,0], # POL
- [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], # ACL
- [2,2,2,2,2,2,2,3,2,2], # LCL
- [0,0,0,0,0,0,0,0,0,0]] # POP
-
-ligaments_fib = [[0,0,0,0,0,0,0,0,0,0], # PCL
- [0,0,0,0,0,0,0,0,0,0], # MCLp
- [0,0,0,0,0,0,0,0,0,0], # MCLd
- [0,0,0,0,0,0,0,0,0,0], # MCLd2
- [0,0,0,0,0,0,0,0,0,0], # POL
- [0,0,0,0,0,0,0,0,0,0], # POL2
- [0,0,0,0,0,0,0,0,0,0], # POL3
- [0,0,0,0,0,0,0,0,0,0], # POL4
- [0,0,0,0,0,0,0,0,0,0], # ACL
- [2,2,2,2,2,2,2,3,2,2], # LCL
- [0,0,0,0,0,0,0,0,0,0]] # POP
-
-for segment in segments:
- if segment == 'femur':
- ligaments = ligaments_fem
- elif segment == 'fibula':
- ligaments = ligaments_fib
+sides = ['L']
+segments = ['femur','tibia', 'fibula']
+short_ssm = [0,1,0]
+no_particles = [4096,2048,4096]
+
+for ind, segment in enumerate(segments):
+ if short_ssm[ind]:
+ short = '_short'
else:
- ligaments = ligaments_tib
+ short = ''
- SSMpoints = [[] for i in range(11)]
- for ind in range(0,11):
- SSMpoints[ind] = [[] for i in range(10)]
+ # Load which SSM points are related to which ligaments
+ 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']]
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'
+ if sides[ind] == 'R':
+ side = '_R'
reflect = ''
else:
- side = 'L'
+ side = '_L'
reflect = '.reflect'
"""SSM part"""
# files from SSM workflow shapeworks
+ particle_file_name = subject + side + '_' + segment + reflect + '.isores.pad.com.center.aligned.clipped.cropped.tpSmoothDT_local.particles'
+ shape_model_folder = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone' + short + r'\new_bone_mri\shape_models'
- file_com = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\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\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\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\new_bone_mri\groomed\com_aligned\S0_' + side + '_' + segment + reflect + '.isores.pad.com.nrrd'
- if segment == 'tibia':
- 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'
- else:
- particle_file = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone\new_bone_mri\shape_models\2048\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\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'
+ path_bones = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output/' + segment + r'_bone' + short + r'\new_bone_mri\input'
mesh_inp = r'C:\Users\mariskawesseli\Documents\GitLab\knee_ssm\OAI\Output\fibula_bone\new_bone_mri\input\S0_L_fibula.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']
-
- # 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)
+ # Create xyz file from particles file
+ pre, ext = os.path.splitext(particle_file_name)
+ particle_file = os.path.join(shape_model_folder, str(no_particles[ind]), particle_file_name)
+ xyz_file = os.path.join(shape_model_folder, pre + '.xyz')
+ copyfile(particle_file, xyz_file)
- # 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(mesh_inp)
- # geometric_measures1 = ms6.apply_filter('compute_geometric_measures')
- # ms6.load_new_mesh(mesh3)
- # geometric_measures2 = ms6.apply_filter('compute_geometric_measures')
- # if reflect == '.reflect':
- # ms6.apply_filter('transform_scale_normalize',axisx=-1,uniformflag=False)
- # translate = geometric_measures1['barycenter']-geometric_measures2['barycenter']
- # 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')
-
- mesh_inp = trimesh.load_mesh(path + '/bone_fibula_1_tissue_rot.stl')
- points_xyz = trimesh.load_mesh(mesh3)
+ # Reflect (mirror) points if needed and translate to the position of the original mesh
+ mesh_inp = trimesh.load_mesh(mesh_inp)
+ points_xyz = trimesh.load_mesh(xyz_file)
if reflect == '.reflect':
M = trimesh.transformations.scale_and_translate((-1, 1, 1))
points_xyz.apply_transform(M)
- Rx = trimesh.transformations.rotation_matrix(90, [0,0,1])
- points_xyz.apply_transform(Rx)
translate = mesh_inp.center_mass-points_xyz.centroid
points_xyz.apply_transform(trimesh.transformations.translation_matrix(translate))
- np.savetxt(path + '\SSM_' + segment + '_transform.xyz', points_xyz.vertices, delimiter=" ")
+ # np.savetxt(path + '\SSM_' + segment + '_transform.xyz', points_xyz.vertices, delimiter=" ") # save intermediate translation to check
# run ICP to get final position SSM point cloud on original mesh
- mesh = trimesh.load_mesh(path + '/bone_fibula_1_tissue_rot.stl') # bone_femur2_2_bone_rot.stl bone_tibia_2_bone_rot
- # 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 + '.npz') # _short
- 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)]
+ icp = trimesh.registration.icp(points_xyz.vertices,mesh_inp,initial=np.identity(4),threshold=1e-5,max_iterations=40,**kwargs)
+ # icp = trimesh.registration.icp(points_xyz.vertices, mesh_inp, initial=icp[0], threshold=1e-5, max_iterations=20,**kwargs) # run icp twice to improve fit
+ points_xyz.apply_transform(icp[0])
+ np.savetxt(path + r'\SSM_' + segment + '_transform_icp.xyz', points_xyz.vertices, delimiter=" ") # save position SSM points on original mesh
+
+ # link which SSM points are related to which ligaments to new point cloud
+ pred_lig_points = points_xyz.vertices[np.hstack(occurances).astype(int)]
+ np.savetxt(path + r'\SSM_' + segment + '_short_pred_points.xyz', np.asarray(pred_lig_points), delimiter=" ")
+ # pred_lig_points_color = np.c_[points_xyz.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])
-
+ print('processing ' + segment + ' done for ' + subject)
\ No newline at end of file
diff --git a/LigamentStudy/scaleOsim.py b/LigamentStudy/scaleOsim.py
index cb79ec4..41081e5 100644
--- a/LigamentStudy/scaleOsim.py
+++ b/LigamentStudy/scaleOsim.py
@@ -2,9 +2,7 @@ import trimesh
import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt
-# import os
-# os.add_dll_directory(r"C:\OpenSim 4.3\bin")
-# import opensim
+
def interpolate_lig_points(lig_points, no_points, plot = 1):
x = lig_points[:, 0]
@@ -160,95 +158,88 @@ osim_model = r'C:\Users\mariskawesseli\Documents\MOBI\data\S0_2_meniscus_lig.osi
# # 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=" ")
+# # 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')
#
-# 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)
+# kwargs = {"scale": False}
+# icp1 = trimesh.registration.icp(mesh2.vertices,mesh1,initial=np.identity(4),threshold=1e-5,max_iterations=20,**kwargs)
#
-# osim_points_tib = 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)],1)
-# np.savetxt(r'C:\Users\mariskawesseli\Documents\LigamentStudy\MRI\S0_prelim\tibia\SSM_tibia_pred_points_osim_interp.xyz', osim_points_tib, delimiter=" ")
+# 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')
#
-# # 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()
+# 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)
+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)
+
+osim_points_tib = 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)],1)
+np.savetxt(r'C:\Users\mariskawesseli\Documents\LigamentStudy\MRI\S0_prelim\tibia\SSM_tibia_pred_points_osim_interp.xyz', osim_points_tib, delimiter=" ")
+
"""fibula"""
# run ICP to get final position SSM point cloud on original mesh
+# mesh OpenSim model - make sure this is high quality
mesh1 = trimesh.load_mesh('C:\opensim-jam\jam-resources\jam-resources-main\models\knee_healthy\smith2019\Geometry\smith2019-R-fibula-bone_remesh.stl')
+# mesh segmented from MRI
mesh2 = trimesh.load_mesh(r'C:\Users\mariskawesseli\Documents\LigamentStudy\MRI\S0_prelim\bone_fibula_1_tissue_rot_remesh.STL')
-# points = trimesh.load_mesh(path + '\SSM_' + segment + '_transform.xyz')
+# Mirror if needed (only for left as SSM/model is right? - check how to deal with left model)
M = trimesh.transformations.scale_and_translate((1,-1,1))
mesh2.apply_transform(M)
+# Rotate segmented bone (check why needed?)
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)
Rz = trimesh.transformations.rotation_matrix(180/(180/np.pi), zaxis)
R = trimesh.transformations.concatenate_matrices(Ry, Rx, Rz)
mesh2.apply_transform(R)
+# Translate segmented mesh to OpenSim bone location
T = trimesh.transformations.translation_matrix(mesh1.center_mass-mesh2.center_mass)
mesh2.apply_transform(T)
mesh2.export(r'C:\Users\mariskawesseli\Documents\LigamentStudy\MRI\S0_prelim\bone_fibula_1_tissue_rot_remesh_test.STL')
+# ICP to fit segmented bone to OpenSim mesh
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_fibula_1_tissue_rot_remesh_icp.STL')
@@ -261,7 +252,6 @@ ligs.apply_transform(T)
ligs.apply_transform(icp[0])
np.savetxt(r'C:\Users\mariskawesseli\Documents\LigamentStudy\MRI\S0_prelim\fibula\SSM_fibula_pred_points_osim.xyz', ligs.vertices, delimiter=" ")
-
ligs = trimesh.load_mesh(r'C:\Users\mariskawesseli\Documents\LigamentStudy\MRI\S0_prelim\fibula\SSM_fibula_pred_points_osim.xyz')
LCL = ligs.vertices[0:79]
@@ -269,10 +259,3 @@ LCL_osim = interpolate_lig_points(LCL,4)
osim_points_tib = np.concatenate([np.asarray(LCL_osim)],1)
np.savetxt(r'C:\Users\mariskawesseli\Documents\LigamentStudy\MRI\S0_prelim\fibula\SSM_fibula_pred_points_osim_interp.xyz', osim_points_tib, 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_osim[0], LCL_osim[1], LCL_osim[2], 'bo')
-# fig2.show()
-# plt.show()
\ No newline at end of file
--
GitLab