Merge pull request #41 from zasexton/main

fixing mis-matched "GlobalNodeID" type expectations

Author: zasexton

svv/simulation/simulation.py

@@ -200,9 +200,9 @@ def build_meshes(self, fluid=True, tissue=False, hausd=0.0001, hsize=None, minra
                     hmin = ((4.0*low_tri_area)/3.0**0.5) ** (0.5)
                     upper_tri_area = area / lower_num_triangles
                     hmax = ((4.0*upper_tri_area)/3.0**0.5) ** (0.5)
-                    #tissue_domain = remesh_surface(tissue_domain, hausd=hausd)
-                #else:
-                #    #tissue_domain = remesh_surface(tissue_domain, hausd=hausd)
+                    tissue_domain = remesh_surface(tissue_domain, hausd=hausd)
+                else:
+                    tissue_domain = remesh_surface(tissue_domain, hausd=hausd)
                 tet_tissue = tetgen.TetGen(tissue_domain)
                 if not fluid:
                     self.synthetic_object.data[0, 0:3] += root_extension * self.synthetic_object.data.get('w_basis', 0)

svv/simulation/utils/extract_faces.py

@@ -407,6 +407,7 @@ def compute_circularity(loop_polydata):
                 npts = face.GetNumberOfPoints()
                 # Canonicalize face nodes to a sorted tuple so orientation/order doesn't matter
                 key = tuple(sorted(face.GetPointId(k) for k in range(npts)))
+                print("key: {} -> {}".format(key, i))
                 face_to_cell[key] = i
     #for i, cap in enumerate(iscap):
     #    if not cap == 1:
@@ -422,9 +423,9 @@ def compute_circularity(loop_polydata):
         if not isinstance(mesh, type(None)):
             wall_surface.point_data["GlobalNodeID"] = numpy.zeros(wall_surface.n_points, dtype=numpy.int32)
             wall_surface.cell_data['GlobalElementID'] = numpy.zeros(wall_surface.n_cells, dtype=numpy.int32)
-            wall_surface.cell_data['ModelFaceID'] = numpy.ones(wall_surface.n_cells, dtype=numpy.int32)
+            wall_surface.cell_data['ModelFaceID'] = numpy.ones(wall_surface.n_cells, dtype=numpy.int32) * i
             _, indices = global_node_tree.query(wall_surface.points)
-            wall_surface.point_data["GlobalNodeID"] = indices.astype(int)
+            wall_surface.point_data["GlobalNodeID"] = indices.astype(numpy.int32)
             # Assign Global Element IDs
             wall_faces = wall_surface.point_data["GlobalNodeID"][wall_surface.faces]
             wall_faces = wall_faces.reshape(-1, 4)[:, 1:].tolist()
@@ -465,10 +466,10 @@ def compute_circularity(loop_polydata):
         if not isinstance(mesh, type(None)):
             cap_surface.point_data["GlobalNodeID"] = numpy.zeros(cap_surface.n_points, dtype=numpy.int32)
             cap_surface.cell_data["GlobalElementID"] = numpy.zeros(cap_surface.n_cells, dtype=numpy.int32)
-            cap_surface.cell_data["ModelFaceID"] = numpy.ones(cap_surface.n_cells, dtype=numpy.int32)
+            cap_surface.cell_data["ModelFaceID"] = numpy.ones(cap_surface.n_cells, dtype=numpy.int32) * (len(walls) + i)
             # Assign Global Node IDs
             _, indices = global_node_tree.query(cap_surface.points)
-            cap_surface.point_data["GlobalNodeID"] = indices.astype(int)
+            cap_surface.point_data["GlobalNodeID"] = indices.astype(numpy.int32)
             # Assign Global Element IDs
             cap_faces = cap_surface.point_data["GlobalNodeID"][cap_surface.faces]
             cap_faces = cap_faces.reshape(-1, 4)[:, 1:].tolist()
@@ -507,10 +508,10 @@ def compute_circularity(loop_polydata):
         if not isinstance(mesh, type(None)):
             lumen_surface.point_data["GlobalNodeID"] = numpy.zeros(lumen_surface.n_points, dtype=int)
             lumen_surface.cell_data["GlobalElementID"] = numpy.zeros(lumen_surface.n_cells, dtype=int)
-            lumen_surface.cell_data["ModelFaceID"] = numpy.ones(lumen_surface.n_cells, dtype=int) * (len(caps) + i + 2)
+            lumen_surface.cell_data["ModelFaceID"] = numpy.ones(lumen_surface.n_cells, dtype=int) * (len(walls) + len(caps) + i)
             # Assign Global Node IDs
             _, indices = global_node_tree.query(lumen_surface.points)
-            lumen_surface.point_data["GlobalNodeID"] = indices.astype(int)
+            lumen_surface.point_data["GlobalNodeID"] = indices.astype(numpy.int32)
             # Assign Global Element IDs
             lumen_faces = lumen_surface.point_data["GlobalNodeID"][lumen_surface.faces]
             lumen_faces = lumen_faces.reshape(-1, 4)[:, 1:].tolist()

svv/tree/branch/bifurcation.py

@@ -1904,6 +1904,19 @@ def cost(x, func=tree_cost_2, d_min=d_min, terminal=terminal,
         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)])
+        vec1 = distal - x
+        vec2 = terminal - x
+        vec3 = proximal - x
+        vec1 = vec1/numpy.linalg.norm(vec1)
+        vec2 = vec2/numpy.linalg.norm(vec2)
+        vec3 = vec3/numpy.linalg.norm(vec3)
+        angle1 = numpy.arctan2(numpy.dot(vec1, vec3))*(180/numpy.pi)
+        angle2 = numpy.arctan2(numpy.dot(vec2, vec3))*(180/numpy.pi)
+        #angle3 = numpy.arccos(numpy.dot(vec3, vec1))*(180/numpy.pi)
+        if angle1 > 90 or angle2 > 90:
+            angle_penalty = penalty
+        else:
+            angle_penalty = 0.0
         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
@@ -1924,7 +1937,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 + triad_penalty, nan=2 * scale + penalty), 0, 2 * scale + penalty)
+        return -1 / np.clip(numpy.nan_to_num(results + triad_penalty + angle_penalty, nan=2 * scale + penalty), 0, 2 * scale + penalty)
         #return results
         #return results
         #return value