Merge pull request #39 from zasexton/main

improving automated 3d mesh generation

Author: zasexton

svv/__init__.py

@@ -1,4 +1,5 @@
-__version__ = "0.0.35"
+__version__ = "0.0.36"
+
 
 # If the optional companion package with compiled accelerators is installed
 # (svv-accelerated), prefer those modules transparently by aliasing them into

svv/simulation/simulation.py

@@ -92,7 +92,10 @@ def build_meshes(self, fluid=True, tissue=False, hausd=0.0001, hsize=None, minra
                             extension_scale += 1.0
                     root_extension = self.synthetic_object.data[0, 21] * extension_scale
                     self.synthetic_object.data[0, 0:3] += root_extension * self.synthetic_object.data.get('w_basis', 0)
+                print("Unioning water-tight model")
                 fluid_surface_mesh = self.synthetic_object.export_solid(watertight=True)
+                print("Finished Unioning water-tight model")
+                print("Tetrahedralizing")
                 tet_fluid = tetgen.TetGen(fluid_surface_mesh)
                 try:
                     tet_fluid.make_manifold(verbose=True)
@@ -182,7 +185,11 @@ def build_meshes(self, fluid=True, tissue=False, hausd=0.0001, hsize=None, minra
                     else:
                         fluid_surface_boolean_mesh = deepcopy(self.fluid_domain_wall_layers[-1])
                 hsize = fluid_surface_boolean_mesh.cell_data["hsize"][0]
-                tissue_domain = remesh_surface(self.synthetic_object.domain.boundary, hausd=hausd) # Check if this should be remeshed
+                try:
+                    tissue_domain = remesh_surface(self.synthetic_object.domain.boundary, hausd=hausd) # Check if this should be remeshed
+                except:
+                    print("REMESHING FAILS: CHECKING FOR TRIANGLE INTERSECTIONS")
+                    tmp_boundary = pymeshfix.MeshFix(self.synthetic_object.domain.boundary)
                 area = tissue_domain.area
                 tissue_domain = boolean(tissue_domain, fluid_surface_boolean_mesh, operation='difference')
                 if fluid:

svv/tree/branch/bifurcation.py

@@ -91,7 +91,7 @@ def add_vessel(tree, **kwargs):
     nonconvex_sampling = kwargs.pop('nonconvex_sampling', 10)
     homogeneous = kwargs.pop('homogeneous', True)
     use_brute = kwargs.pop('use_brute', False)
-    max_iter = kwargs.pop('max_iter', 10)
+    max_iter = kwargs.pop('max_iter', 100)
     return_cost = kwargs.pop('return_cost', False)
     #defualt_threshold = ((tree.domain.mesh.volume ** (1/3)) /
     #                     (tree.n_terminals ** threshold_exponent)) + tree.data[0, 21]*2.0
@@ -218,8 +218,8 @@ def callback(xk):
                                 #plt.colorbar(label='Function values')
                                 #plt.show()
                                 cons = [{"type": "ineq", "fun": lambda a: 1 - a[0] - a[1]}]
-                                result = minimize(cost, x0, bounds=[(0, 1.0), (0, 1.0)], callback=callback,
-                                                  options={'maxiter':max_iter})
+                                result = minimize(cost, x0, bounds=[(0.05, 0.95), (0.05, 0.95)], callback=callback,
+                                                  options={'maxiter':max_iter},constraints=cons, method="COBYLA")
                                 #print('SOLUTION: {}'.format(result.x))
                                 #print('SOLUTION FUN: {}'.format(result.fun))
                                 bifurcation_point = triad(result.x)
@@ -1901,10 +1901,10 @@ def cost(x, func=tree_cost_2, d_min=d_min, terminal=terminal,
              scale=tree_scale,connectivity=tree.connectivity):
         x = triad(x)
         #dists = close_exact_point(lines, x)
-        #dists = numpy.array([numpy.linalg.norm(lines[0, 0:3] - x),
-        #                     numpy.linalg.norm(lines[0, 3:6] - x),
-        #                     numpy.linalg.norm(lines[1, 3:6] - x)])
-        #triad_penalty = numpy.max([0.0, -1.0 * numpy.min(dists - d_min)])/d_min * penalty
+        dists = numpy.array([numpy.linalg.norm(lines[0, 0:3] - x),
+                             numpy.linalg.norm(lines[0, 3:6] - x),
+                             numpy.linalg.norm(lines[1, 3:6] - x)])
+        triad_penalty = numpy.max([0.0, -1.0 * numpy.min(dists - d_min)])/d_min * penalty
         #[TODO] angle penalty
         #[TODO] require that resulting parent vessel is at least a certain length? remove buffer region around triad
         # points
@@ -1924,7 +1924,7 @@ def cost(x, func=tree_cost_2, d_min=d_min, terminal=terminal,
         #assert results > tree_scale, '{} results < {} tree_scale'.format(results, tree_scale)
         #return (((np.clip(numpy.nan_to_num(results - scale, nan=2*scale+penalty), 0, 2*scale+penalty) + triad_penalty))/(scale+penalty))# + 1.0
         #return -1/np.clip(numpy.nan_to_num(results - scale, nan=2*scale+penalty), 0, 2*scale+penalty)
-        return -1 / np.clip(numpy.nan_to_num(results, nan=2 * scale + penalty), 0, 2 * scale + penalty)
+        return -1 / np.clip(numpy.nan_to_num(results + triad_penalty, nan=2 * scale + penalty), 0, 2 * scale + penalty)
         #return results
         #return results
         #return value