import sys import os import inspect import time import subprocess import glob try: from lxml import etree as et # not in my standard Salome's python but better except ImportError: import xml.etree.cElementTree as et #import numpy as np try: import salome import salome_notebook import SMESH, SALOMEDS from salome.smesh import smeshBuilder except ImportError: print('error in imports') sys.exit() try: from MEDLoader import MEDLoader except ImportError: import MEDLoader def clearSalome(nb): '''Delete all meshes in study''' for compName in ["SMESH", "GEOM"]: comp = salome.myStudy.FindComponent(compName) if comp: iterator = salome.myStudy.NewChildIterator(comp) while iterator.More(): sobj = iterator.Value() iterator.Next() nb.RemoveObjectWithChildren(sobj) def closeSalome(): """Try to close out of Salome gracefully""" # there is nothing graceful about this try: sys.stdout = 'redirect to nowhere...' # so I don't see the junk it prints from killSalomeWithPort import killMyPort killMyPort(os.getenv('NSPORT')) except: pass def modelSeparation(smesh, ams_path): # Make a new directory to hold the sliced mesh files and various outputs split_dir = ams_path.rsplit(os.path.sep, 1)[0] + os.path.sep print(split_dir) new_dir = split_dir + 'SlicedMeshes'+ os.path.sep print(new_dir) if os.path.isdir(new_dir): pass else: os.mkdir(new_dir) #os.chdir(new_dir) files = glob.glob(new_dir + "*") for f in files: os.remove(f) # Export the 3 volumes as separate meshes ([Flesh], status) = smesh.CreateMeshesFromMED(ams_path) #[Bone, Flesh_All_Volumes, Muscle, Fat, Flesh_All_Faces, Flesh__Probe_ContactFaces] = Flesh.GetGroups() [Flesh__Probe_ContactFaces, Flesh_All_Faces, Fat, Flesh_All_Volumes, Muscle, Bone] = Flesh.GetGroups() try: print(new_dir + 'BoneLinear.med') # Flesh.ExportMED(new_dir + 'BoneLinear.med', 0, 40, 1, Bone, 1, [], '', -1) # Flesh.ExportMED(new_dir + 'FatLinear.med', 0, 40, 1, Fat, 1, [], '', -1) # Flesh.ExportMED(new_dir + 'MuscleLinear.med', 0, 40, 1, Muscle, 1, [], '', -1) Flesh.ExportMED(new_dir + 'BoneLinear.med', 0, SMESH.MED_V2_2, 1, Bone, 1,[], '') Flesh.ExportMED(new_dir + 'MuscleLinear.med', 0, SMESH.MED_V2_2, 1, Muscle, 1, [], '') Flesh.ExportMED(new_dir + 'FatLinear.med', 0, SMESH.MED_V2_2, 1, Fat, 1, [], '') except: print('ExportPartToMED() failed. Invalid file name?') return new_dir def meshCut(salome_path, start_med, end_med, mesh_name, top_cut, bottom_cut, plane_norm, plane_vertex): '''Call the smesh plugin, MeshCut, from the terminal since it has no salome command to plane cut a mesh. Convert all elements to tets after cutting. Would we ever want to change the tolerance? Terminal command takes the form: MeshCut BoneLinear.med Output.med MeshName Above Below 0 0 -1 0, 0 -930 0.1''' #FOR WINDOWS #cut_command = salome_path + os.path.join("MODULES","SMESH","RELEASE","SMESH_INSTALL","bin","salome","MeshCut.exe") # FOR LINUX: r'/home/doherts/salome/appli_V9_3_0/bin/salome/MeshCut' #FOR LINUX cut_command = salome_path + os.path.join("bin", "salome", "MeshCut") #print(cut_command) if os.path.isfile(cut_command): meshcut_str = (13*"{} ").format(cut_command, start_med, end_med, mesh_name, top_cut,bottom_cut, plane_norm[0], plane_norm[1], plane_norm[2], plane_vertex[0], plane_vertex[1], plane_vertex[2],.1) proc = subprocess.Popen(meshcut_str, shell=True, stdout=subprocess.PIPE) print(proc.communicate()) else: print("Is your meshcut salome path correct? Check meshCut function and the cut_command variable") #Checking for MODULES\\INSTALL\\SMESH\\RELEASE\\SMESH_INSTALL\\bin\\salome\\MeshCut.exe") #/home/doherts/salome/appli_V9_3_0/bin/salome/MeshCut") def tetClean(nb, smesh, med_file): '''Remove extraneous groups and ensure the mesh is entirely tetrahedrals''' ([CutMesh], status) = smesh.CreateMeshesFromMED(med_file) group_names = CutMesh.GetGroupNames() print(group_names) # print GroupNames[0].replace(" ", "") if len(group_names) == 2: group_one = group_names[0].replace(" ", "") print (group_names[0].replace(" ", "") == "keep") if group_one == 'keep': try: [keep, delete] = CutMesh.GetGroups() CutMesh.SplitVolumesIntoTetra(CutMesh, 1) #CutMesh.ExportMED(med_file, 0, 40, 1, keep, 1, [], '', -1) #os.remove(med_file) CutMesh.ExportMED(med_file, 0, SMESH.MED_V2_2, 1, keep, 1, [], '') #CutMesh.ExportMED(med_file, 0, SMESH.MED_V2_2, 1, delete, 1, [], '') except: pass else: try: [delete, keep] = CutMesh.GetGroups() CutMesh.SplitVolumesIntoTetra(CutMesh, 1) #CutMesh.ExportMED(med_file, 0, 40, 1, keep, 1, [], '', -1) #os.remove(med_file) CutMesh.ExportMED(med_file, 0, SMESH.MED_V2_2, 1, keep, 1, [], '') except: pass else: print("No cut happened") pass clearSalome(nb) def meshCutFolder(salome_path, nb, smesh, folder, contains_name, norm, vertex): '''Cut all meshes in a folder with the same plane. Control which meshes to cut by the name variable. ''' for file in glob.glob(folder + '*'): print(file) if '.med' in file and contains_name in file: if 'cut' not in file: print("starting meshcut") new_file = file.split('.')[0] + '_cut.med' meshCut(salome_path,file, new_file, 'CutMesh', 'keep', 'delete', norm, vertex) tetClean(nb, smesh, new_file) else: new_file = file.split('.')[0] + '_temp.med' meshCut(salome_path, file, new_file, 'CutMesh', 'keep', 'delete', norm, vertex) tetClean(nb, smesh, new_file) #os.remove(file) os.rename(new_file, new_file.split('_temp.med')[0] + '.med') else: pass def findProbeMotion(febioFile): '''Read and return the vector that prescribes the probes displacement''' f = open(febioFile, 'r') tree = et.parse(f) febio_root = tree.getroot() BCs = febio_root.find('Boundary').find('rigid_body').findall('prescribed') xFebDisp = float(BCs[0].text) yFebDisp = float(BCs[1].text) zFebDisp = float(BCs[2].text) probeVector = [xFebDisp,yFebDisp,zFebDisp] f.close() return probeVector def probeMotionModel(nb, smesh, region, salome_path, new_dir, probe_mesh, d_scale, t_scale, s_scale): '''This really just slightly shifts the location of the cropped model to account for probe motion''' #13384 was previously found to be the center of the probes indentation face (salome is 1 indexed) center_node = 13384 x, y, z = probe_mesh.GetNodeXYZ(center_node) center_faces = probe_mesh.GetNodeInverseElements(center_node) center_norms = [] for face in center_faces: center_norms.append(probe_mesh.GetFaceNormal(face, normalized=True)) print("Probe Normal Axis: " , center_norms) xnorm = [] ynorm = [] znorm = [] for norm in center_norms: xnorm.append(norm[0]) ynorm.append(norm[1]) znorm.append(norm[2]) average_norm = [sum(xnorm) / len(xnorm), sum(ynorm) / len(ynorm), sum(znorm) / len(znorm)] print("Probe Normal Axis: ", average_norm) clearSalome(nb) # these values were generated from other code, as a linear fit of the registration data scaling_factor = 2 # THIS IS THE DEFAULT I USED, CHANGE IF YOU USE A DIFFERENT PROBE MOTION region_dict = {} region_dict['008UA'] = [2.0690775274597475, -4.680213714438663, 0.33758168068581673] region_dict['008LA'] = [-0.37713671420684397, -2.947740900081625, -0.4728674816914108] region_dict['008UL'] = [3.8011597697464152, -5.0758354202077305, 0.5264698493663289] region_dict['008LL'] = [-0.8741175625230704, -2.480189901433061, 1.3337193437535921] region_dict['006UA'] = [-2.860630733468474, -7.698406183503814, 2.119024510460199] region_dict['006LA'] = [1.3591328966912897, -5.727680728681045, -0.20363212279028708] region_dict['006UL'] = [-2.6967931657448583, -2.0569732451906075, 0.9117781511680523] region_dict['006LL'] = [-0.08612317828278482, -4.396007117918238, 0.47309440274900805] vector = region_dict[region] x = x + vector[0]*scaling_factor # NO CHANGE TO Y AS THAT IS THE DEPTH z = z + vector[2]*scaling_factor #export dynamic library paths to enable running meshcut from the command line if os.path.isdir(salome_path): #Assume the mesh runs along the z-axis meshCutFolder(salome_path, nb, smesh, new_dir, 'Linear', [0, 0, 1], [0, 0, z-10*t_scale]) meshCutFolder(salome_path, nb, smesh, new_dir, 'cut', [1, 0, 0], [x-10*s_scale, 0, 0]) meshCutFolder(salome_path, nb, smesh, new_dir, 'cut', [-1, 0, 0], [x+10*s_scale, 0, 0]) meshCutFolder(salome_path, nb, smesh, new_dir, 'cut', [0, 0, -1], [0, 0, z+10*t_scale]) meshCutFolder(salome_path, nb, smesh, new_dir, 'cut', [-average_norm[0], -average_norm[1], -average_norm[2]], [x + 10*d_scale*average_norm[0], y + 10*d_scale*average_norm[1], z + 10*d_scale*average_norm[2]]) #meshCutFolder(nb, smesh, new_dir, 'cut', []) clearSalome(nb) else: print("Error when using MeshCutFolder") #/bin/salome/MeshCut ([Bone], status) = smesh.CreateMeshesFromMED(new_dir + 'BoneLinear_cut.med') ([Fat], status) = smesh.CreateMeshesFromMED(new_dir + 'FatLinear_cut.med') ([Muscle], status) = smesh.CreateMeshesFromMED(new_dir + 'MuscleLinear_cut.med') Compound_Mesh_1 = smesh.Concatenate([Bone.GetMesh(), Fat.GetMesh(), Muscle.GetMesh()], 1, 1, 1e-005, True) [Bone_Nodes, Bone_Volumes, Fat_Nodes, Fat_Volumes, Muscle_Nodes, Muscle_Volumes] = Compound_Mesh_1.GetGroups() #Set names of Mesh objects smesh.SetName(Compound_Mesh_1.GetMesh(), 'Cropped_Mesh') smesh.SetName(Bone_Nodes, 'Bone_Nodes') smesh.SetName(Fat_Nodes, 'Fat_Nodes') smesh.SetName(Muscle_Nodes, 'Muscle_Nodes') smesh.SetName(Bone_Volumes, 'Bone') smesh.SetName(Fat_Volumes, 'Fat') smesh.SetName(Muscle_Volumes, 'Muscle') Compound_Mesh_1.ExportMED(new_dir + 'Cropped_Model.med', 0, SMESH.MED_V2_2, 1, None, 1) print('\n\nFiles will be located at {}\n\n'.format(new_dir)) def probeLineOfActionModel(nb, smesh, salome_path, new_dir, probe_mesh, d_scale, t_scale, s_scale): #13384 was previously found to be the center of the probes indentation face (salome is 1 indexed) center_node = 13384 x, y, z = probe_mesh.GetNodeXYZ(center_node) center_faces = probe_mesh.GetNodeInverseElements(center_node) center_norms = [] for face in center_faces: center_norms.append(probe_mesh.GetFaceNormal(face, normalized=True)) print("Probe Normal Axis: " , center_norms) xnorm = [] ynorm = [] znorm = [] for norm in center_norms: xnorm.append(norm[0]) ynorm.append(norm[1]) znorm.append(norm[2]) average_norm = [sum(xnorm) / len(xnorm), sum(ynorm) / len(ynorm), sum(znorm) / len(znorm)] print("Probe Normal Axis: ", average_norm) clearSalome(nb) #export dynamic library paths to enable running meshcut from the command line if os.path.isdir(salome_path): #Assume the mesh runs along the z-axis #FIX: Can i optimize plane norms? meshCutFolder(salome_path, nb, smesh, new_dir, 'Linear', [0, 0, 1], [0, 0, z-10*t_scale]) meshCutFolder(salome_path, nb, smesh, new_dir, 'cut', [1, 0, 0], [x-10*s_scale, 0, 0]) meshCutFolder(salome_path, nb, smesh, new_dir, 'cut', [-1, 0, 0], [x+10*s_scale, 0, 0]) meshCutFolder(salome_path, nb, smesh, new_dir, 'cut', [0, 0, -1], [0, 0, z+10*t_scale]) meshCutFolder(salome_path, nb, smesh, new_dir, 'cut', [-average_norm[0], -average_norm[1], -average_norm[2]], [x + 10*d_scale*average_norm[0], y + 10*d_scale*average_norm[1], z + 10*d_scale*average_norm[2]]) #meshCutFolder(nb, smesh, new_dir, 'cut', []) clearSalome(nb) else: print("Error with salome path") #/bin/salome/MeshCut ([Bone], status) = smesh.CreateMeshesFromMED(new_dir + 'BoneLinear_cut.med') ([Fat], status) = smesh.CreateMeshesFromMED(new_dir + 'FatLinear_cut.med') ([Muscle], status) = smesh.CreateMeshesFromMED(new_dir + 'MuscleLinear_cut.med') Compound_Mesh_1 = smesh.Concatenate([Bone.GetMesh(), Fat.GetMesh(), Muscle.GetMesh()], 1, 1, 1e-005, True) [Bone_Nodes, Bone_Volumes, Fat_Nodes, Fat_Volumes, Muscle_Nodes, Muscle_Volumes] = Compound_Mesh_1.GetGroups() #Set names of Mesh objects smesh.SetName(Compound_Mesh_1.GetMesh(), 'Cropped_Mesh') smesh.SetName(Bone_Nodes, 'Bone_Nodes') smesh.SetName(Fat_Nodes, 'Fat_Nodes') smesh.SetName(Muscle_Nodes, 'Muscle_Nodes') smesh.SetName(Bone_Volumes, 'Bone') smesh.SetName(Fat_Volumes, 'Fat') smesh.SetName(Muscle_Volumes, 'Muscle') Compound_Mesh_1.ExportMED(new_dir + 'Cropped_Model.med', 0, SMESH.MED_V2_2, 1, None, 1) print('\n\nFiles will be located at {}\n\n'.format(new_dir)) def layeredCreation(nb, smesh, region, ams_path, salome_path, probe_path, d_scale, t_scale, s_scale, probe_direction_bool): '''Create a layered cube model''' #Separate the adaptive meshing model into 3 components (fat, bone, muscle) new_dir = modelSeparation(smesh, ams_path) #Find the probe bounding box #print(path+'{}_ProbeFiles\\9L4_transformed_{}_rmOverlap2.stl'.format(multis_number,multis_region)) probe_mesh = smesh.CreateMeshesFromSTL(probe_path) # bounding_box = probe_mesh.BoundingBox() # print("Current Bounding Box: " ,bounding_box) if probe_direction_bool is True: #findProbeMotion() probeMotionModel(nb, smesh, region, salome_path, new_dir, probe_mesh, d_scale, t_scale, s_scale) else: probeLineOfActionModel(nb, smesh, salome_path, new_dir, probe_mesh,d_scale, t_scale, s_scale) if __name__ == '__main__': # multis_number = '008' # multis_region = 'LowerLeg' # region = multis_number+ ''.join([c for c in multis_region if c.isupper()]) # # path = os.path.join("C:/","Users","sbd","Desktop","WorkStuff","Documents","MULTIS") + os.path.sep # probe_path = path+ os.path.join("008_ProbeFiles", "9L4_transformed_LowerLeg_rmOverlap2.stl") # change this to probe path "{}_ProbeFiles", "9L4_transformed_{}_rmOverlap2.stl" # adaptive_mesh_path = path + os.path.join("008_LL","CMULTIS008-1_LL_adapt.med") # change this to adaptive mesh path # salome_path = os.path.join("C:/","Users","sbd","SALOME-8.3.0-WIN64") + os.path.sep multis_number = '006' multis_region = 'UpperLeg' med_file = 'leg_3.med' multis_abbreviation = ''.join([c for c in multis_region if c.isupper()]) #i.e. "UL" for UpperLeg region = multis_number+multis_abbreviation path = os.path.join("/home", "doherts", "Documents", "MULTIS") + os.path.sep probe_path = path + os.path.join("{}_ProbeFiles".format(multis_number),"9L4_transformed_{}_AvgPos_rmOverlap2.stl".format(multis_region)) # change this to probe path "{}_ProbeFiles", "9L4_transformed_{}_rmOverlap2.stl" adaptive_mesh_path = path + os.path.join("CroppedLayered", "{}".format(region), med_file) # change this to adaptive mesh path salome_path = os.path.join("/home", "doherts", "salome", "appli_V9_3_0") + os.path.sep depth_scale = 5 #how deep the cut should be along the probes line of action , 5 seems to work reasonably well transverse_scale = 6 # needs to be over ~3 to fit large dimension for face of probe saggital_scale = 3 # needs to be over ~2 to fit smaller dimension for face of probe salome.salome_init() theStudy = salome.myStudy notebook = salome_notebook.NoteBook(theStudy) smesh = smeshBuilder.New(theStudy) nb = theStudy.NewBuilder() layeredCreation(nb, smesh, region, adaptive_mesh_path, salome_path, probe_path, depth_scale, transverse_scale, saggital_scale, probe_direction_bool=True)