1. With the python API while I can generate the model I cannot successfully mesh the model. I get an error about self-intersecting facets (see full print out below). However, I can import the same ctgr files into the simvascular GUI and generate a model there that I can successfully mesh. I'm wondering what my code is missing to replicate the GUI process?
2. Secondly, when I generate the model (Create Model button in GUI) in simvascular it recognizes each individual ctgr group as a wall and gives unique outlet ids. In my python code when I create the model there is only one wall and all the outlets share the same id. How can I replicate the GUI behaviour using the python API?
PYTHON CODE
Code: Select all
import vtk
import os
import sv
def get_profile_contour(contours, cid, npts):
cont = contours[cid]
cont_pd = cont.get_polydata()
cont_ipd = sv.geometry.interpolate_closed_curve(polydata=cont_pd, number_of_points=npts)
return cont_ipd
def loft_contours_polydata(contours):
"""Lofts contours for use in polydata reconstruction"""
num_samples = 50
use_distance = True
aligned_contours = []
for index in range(0, len(contours)):
cont = get_profile_contour(contours, index, num_samples)
if index == 0:
cont_align = cont
else:
cont_align = sv.geometry.align_profile(last_cont_align, cont, use_distance)
last_cont_align = cont_align
aligned_contours.append(cont_align)
options = sv.geometry.LoftOptions()
options.interpolate_spline_points = True
options.num_spline_points = 400 # contorls number of elements along vessel length
options.num_long_points = 200 # controls how well matches initial geometry (more curved vessels requires greater values)
lofted_surface = sv.geometry.loft(polydata_list=aligned_contours, loft_options=options)
# Create a new solid from the lofted solid.
lofted_model = sv.modeling.PolyData()
lofted_model.set_surface(surface=lofted_surface)
return lofted_model
def merge(main_contours, branch_contours, recon_method):
"""Merges branches to main vessel"""
if recon_method == 'vtp':
lofted_contours = loft_contours_polydata(main_contours)
lofted_capped = sv.vmtk.cap(surface=lofted_contours.get_polydata(),
use_center=False)
lofted_model = sv.modeling.PolyData()
lofted_model.set_surface(lofted_capped)
for branch_contour in branch_contours:
# merge vessels polydata
if recon_method == 'vtp':
modeler = sv.modeling.Modeler(sv.modeling.Kernel.POLYDATA)
lofted_contours = loft_contours_polydata(branch_contour)
lofted_capped = sv.vmtk.cap(surface=lofted_contours.get_polydata(),
use_center=False)
lofted_model1 = sv.modeling.PolyData()
lofted_model1.set_surface(lofted_capped)
lofted_model = modeler.union(model1=lofted_model, model2=lofted_model1)
return lofted_model
def read_contours(segmentation_path):
'''Read an SV contour group file.
'''
contour_group = sv.segmentation.Series(segmentation_path)
num_conts = contour_group.get_num_segmentations()
contours = []
for i in range(num_conts):
cont = contour_group.get_segmentation(i)
contours.append(cont)
return contours
def create_mesh(model):
mesher = sv.meshing.TetGen()
caps = model.identify_caps()
wall_id = [i + 1 for i in range(len(caps)) if caps[i]]
print("Wall ids are {}".format(wall_id))
# Set the model for the mesher.
mesher.set_model(model)
mesher.set_walls(wall_id)
face_ids = mesher.get_model_face_ids()
print(face_ids)
edge_size = 1.5
print("Set meshing options ... ")
options = sv.meshing.TetGenOptions(global_edge_size=edge_size, surface_mesh_flag=True, volume_mesh_flag=True)
options.radius_meshing_on = False
options.boundary_layer_inside = True
options.radius_meshing_scale = 0.4
# Generate the mesh.
mesher.generate_mesh(options)
# Get the mesh as a vtkUnstructuredGrid.
mesh = mesher.get_mesh()
print("Mesh:")
print(" Number of nodes: {0:d}".format(mesh.GetNumberOfPoints()))
print(" Number of elements: {0:d}".format(mesh.GetNumberOfCells()))
# Write the mesh.
mesh_series = sv.meshing.Series()
mesh_series.write('mesh.msh')
mesh_path = 'mesh.vtu'
mesher.write_mesh(mesh_path)
if __name__ == "__main__":
os.chdir('Segmentations')
main_contours = read_contours('CX.ctgr')
branch_contours = []
branch_contours.append(read_contours('OM1_scaled.ctgr'))
branch_contours.append(read_contours('OM2_scaled.ctgr'))
recon_method= 'vtp'
model = merge(main_contours, branch_contours, recon_method)
create_mesh(model)
Checking surface mesh
Regions: 1
Number of Free Edges on Surface: 209
Number of Non-Manifold Edges on Surface: 0
Iteration 1/10
Iteration 2/10
Iteration 3/10
Iteration 4/10
Iteration 5/10
Iteration 6/10
Iteration 7/10
Iteration 8/10
Iteration 9/10
Iteration 10/10
Final mesh improvement
Converting to TetGen...
Adding Facet Markers...
TetGen Meshing Started...
Delaunizing vertices...
Delaunay seconds: 0.027234
Creating surface mesh ...
Found two nearly self-intersecting facets.
1st: [41, 39, 40] #2
2nd: [41, 39, 6052] #2
The dihedral angle between them is 0.0150374 degree.
Hint: You may use -p/# to decrease the dihedral angle tolerance 0.1 (degree).
ERROR: TetGen quit and returned error code 3