Added provide_inlet method

morphman/automated_landmarking/automated_landmarking_piccinelli.py

@@ -55,9 +55,9 @@ def landmarking_piccinelli(centerline, base_path, approximation_method, algorith
 
     elif approximation_method == "vmtk":
         line = centerline
-        line_curv = vmtk_compute_geometric_features(centerline, True, outputsmoothed=False,
+        line_curv = vmtk_compute_geometric_features(centerline, True, output_smoothed=False,
                                                     factor=smoothing_factor_curv, iterations=iterations)
-        line_tor = vmtk_compute_geometric_features(centerline, True, outputsmoothed=False,
+        line_tor = vmtk_compute_geometric_features(centerline, True, output_smoothed=False,
                                                    factor=smoothing_factor_torsion, iterations=iterations)
         # Get curvature and torsion, find peaks
         curvature = get_point_data_array("Curvature", line_curv)

morphman/common/centerline_operations.py

@@ -36,7 +36,7 @@ def get_bifurcating_and_diverging_point_data(centerline, centerline_bif, tol):
     cl1 = extract_single_line(centerline, 0)
     cl2 = extract_single_line(centerline, 1)
 
-    # Declear dictionary to hold results
+    # Declare dictionary to hold results
     data = {"bif": {}, 0: {}, 1: {}}
 
     # Find lower clipping point
@@ -117,26 +117,26 @@ def get_manipulated_centerlines(patch_cl, dx, p1, p2, diverging_id, diverging_ce
         locator = get_vtk_point_locator(line)
         id1 = locator.FindClosestPoint(p1)
         if diverging_id is not None and i == (number_of_cells - 1):
-            # Note: Reuse id2 and idmid from previous loop
+            # Note: Reuse id2 and id_mid from previous loop
             pass
         else:
             id2 = locator.FindClosestPoint(p2)
-            idmid = int((id1 + id2) * 0.5)
+            id_mid = int((id1 + id2) * 0.5)
 
         for p in range(line.GetNumberOfPoints()):
             point = line.GetPoint(p)
             cl_id = locator.FindClosestPoint(point)
 
             if direction == "horizont":
                 if diverging_id is not None and i == (number_of_cells - 1) and diverging_id < cl_id:
-                    dist = -dx * (diverging_id - idmid) ** 0.5 / (diverging_id - idmid) ** 0.5
+                    dist = -dx * (diverging_id - id_mid) ** 0.5 / (diverging_id - id_mid) ** 0.5
                 else:
                     if cl_id < id1:
                         dist = dx
-                    elif id1 <= cl_id < idmid:
-                        dist = dx * (idmid ** 2 - cl_id ** 2) / (idmid ** 2 - id1 ** 2)
-                    elif idmid <= cl_id < (id2 - 1):
-                        dist = -dx * (cl_id - idmid) ** 0.5 / (id2 - idmid) ** 0.5
+                    elif id1 <= cl_id < id_mid:
+                        dist = dx * (id_mid ** 2 - cl_id ** 2) / (id_mid ** 2 - id1 ** 2)
+                    elif id_mid <= cl_id < (id2 - 1):
+                        dist = -dx * (cl_id - id_mid) ** 0.5 / (id2 - id_mid) ** 0.5
                     else:
                         dist = -dx
 
@@ -222,15 +222,15 @@ def get_curvilinear_coordinate(line):
         line (vtkPolyData): Input centerline
 
     Returns:
-        curv_coor (ndarray): Array of abscissa points.
+        curvilinear_coordinate (ndarray): Array of abscissa points.
     """
-    curv_coor = np.zeros(line.GetNumberOfPoints())
+    curvilinear_coordinate = np.zeros(line.GetNumberOfPoints())
     for i in range(line.GetNumberOfPoints() - 1):
         pnt1 = np.asarray(line.GetPoints().GetPoint(i))
         pnt2 = np.asarray(line.GetPoints().GetPoint(i + 1))
-        curv_coor[i + 1] = np.sqrt(np.sum((pnt1 - pnt2) ** 2)) + curv_coor[i]
+        curvilinear_coordinate[i + 1] = np.sqrt(np.sum((pnt1 - pnt2) ** 2)) + curvilinear_coordinate[i]
 
-    return curv_coor
+    return curvilinear_coordinate
 
 
 def get_centerline_tolerance(centerline, n=50):
@@ -255,15 +255,15 @@ def get_centerline_tolerance(centerline, n=50):
 
 def get_clipped_diverging_centerline(centerline, clip_start_point, clip_end_id):
     """
-    Clip the opthamlic artery if present.
+    Clip the ophthalmic artery if present.
 
     Args:
-        centerline (vtkPolyData): Line representing the opthalmic artery centerline.
-        clip_start_point (tuple): Point at entrance of opthalmic artery.
-        clip_end_id (int): ID of point at end of opthalmic artery.
+        centerline (vtkPolyData): Line representing the ophthalmic artery centerline.
+        clip_start_point (tuple): Point at entrance of ophthalmic artery.
+        clip_end_id (int): ID of point at end of ophthalmic artery.
 
     Returns:
-        patch_eye (vtkPolyData): Voronoi diagram representing opthalmic artery.
+        patch_eye (vtkPolyData): Voronoi diagram representing ophthalmic artery.
     """
     points = [clip_start_point, centerline.GetPoint(clip_end_id)]
     div_points = vtk.vtkPoints()
@@ -286,7 +286,7 @@ def get_line_to_change(surface, centerline, region_of_interest, method, region_p
         centerline (vtkPolyData): Centerline in geometry.
         region_of_interest (str): Method for setting the region of interest ['manual' | 'commandline' | 'first_line']
         method (str): Determines which kind of manipulation is performed.
-        region_points (list): If region_of_interest is 'commandline', this a flatten list of the start and endpoint.
+        region_points (list): If region_of_interest is 'commandline', a flattened list of the start and endpoint.
         stenosis_length (float): Multiplier used to determine the length of the stenosis-affected area.
 
     Returns:
@@ -488,7 +488,7 @@ def get_region_of_interest_and_diverging_centerlines(centerlines_complete, regio
         region_points (ndarray): Two points determining the region of interest.
 
     Returns:
-        centerlines (vtkPolyData): Centerlines excluding divering lines.
+        centerlines (vtkPolyData): Centerlines excluding diverging lines.
         diverging_centerlines (vtkPolyData): Centerlines diverging from the region of interest.
         region_points (ndarray): Sorted region points.
         region_points_vtk (vtkPoints): Sorted region points as vtkData.
@@ -548,7 +548,7 @@ def compute_discrete_derivatives(line, neigh=10):
 
     Args:
         line (vtkPolyData): Line to compute geometry from.
-        neigh (int): Number of naboring points.
+        neigh (int): Number of neighboring points.
 
     Returns:
         line (vtkPolyData): Output line with geometrical parameters.
@@ -689,7 +689,7 @@ def get_k1k2_basis(curvature, line):
         E1[i, :] = V[:, 1]
         E2[i, :] = V[:, 2]
 
-    # Compute k_1, k_2 furfilling T' = curv(s)*N(s) = k_1(s)*E_1(s) + k_2(s)*E_2(s).
+    # Compute k_1, k_2 fulfilling T' = curv(s)*N(s) = k_1(s)*E_1(s) + k_2(s)*E_2(s).
     # This is simply a change of basis for the curvature vector N. The room
     # of k_1 and k_2 can be used to express both the curvature and the
     # torsion.
@@ -717,11 +717,11 @@ def compute_splined_centerline(line, get_curv=False, isline=False, nknots=50, ge
         get_curv (bool): Computes curvature profile if True.
         isline (bool): Determines if centerline object is a line or points.
         nknots (int): Number of knots.
-        get_stats (bool): Determines if curve attribuites are computed or not.
+        get_stats (bool): Determines if curve attributes are computed or not.
         get_misr (bool): Determines if MISR values are computed or not.
 
     Returns:
-        line (vtkPolyData): Splined centerline data.
+        line (vtkPolyData): Spline of centerline data.
     Returns:
         curv (ndarray): Curvature profile.
     """

morphman/common/common.py

@@ -82,7 +82,7 @@ def get_parameters(folder):
     Returns:
         data (dict): The data in the info file.
     """
-    # If info.json does not exists, return an empty dict
+    # If info.json does not exist, return an empty dict
     if not path.isfile(folder + "_info.json"):
         return {}
 
@@ -121,7 +121,7 @@ def convert_numpy_data_to_polydata(data, header, TNB=None, PT=None):
         PT (numpy.ndarray): Data array.
 
     Returns:
-        line (vtkPolyData): Line couple with all the new data.
+        line (vtkPolyData): Line-couple with all the new data.
     """
     line = vtk.vtkPolyData()
     cell_array = vtk.vtkCellArray()
@@ -504,15 +504,15 @@ def get_rotation_matrix(u, angle):
         angle (float): Angle to rotate in radians
 
     Returns:
-        R (ndarray): Rotation matrix
+        rotation_matrix (ndarray): Rotation matrix
     """
     u_cross_matrix = np.asarray([[0, -u[2], u[1]],
                                  [u[2], 0, -u[0]],
                                  [-u[1], u[0], 0]])
     u_outer = np.outer(u, u)
-    R = np.cos(angle) * np.eye(3) + np.sin(angle) * u_cross_matrix + (1 - np.cos(angle)) * u_outer
+    rotation_matrix = np.cos(angle) * np.eye(3) + np.sin(angle) * u_cross_matrix + (1 - np.cos(angle)) * u_outer
 
-    return R
+    return rotation_matrix
 
 
 def get_angle(a, b):

morphman/common/surface_operations.py

@@ -42,14 +42,15 @@ def get_relevant_outlets(surface, base_path):
     return relevant_outlets
 
 
-def compute_centers(polydata, case_path=None):
+def compute_centers(polydata, case_path=None, select_inlet=False):
     """
     Compute the center of all the openings in the surface. The inlet is chosen based on
     the largest area.
 
     Args:
-        polydata (vtkPolyData): centers of the openings
-        case_path (str): path to case directory.
+        polydata (vtkPolyData): Centers of the openings
+        case_path (str): Path to case directory.
+        select_inlet (bool): Let user select inlet manually
 
     Returns:
         inlet (list): A list of points.
@@ -90,12 +91,18 @@ def compute_centers(polydata, case_path=None):
         delaunay.SetInputData(boundary_points)
         delaunay.Update()
 
-        # Add quanteties
+        # Add quantities
         area.append(vtk_compute_mass_properties(delaunay.GetOutput()))
         center.append(np.mean(tmp_points, axis=0))
 
+    # Select inlet manually or based on area
+    if select_inlet:
+        inlet_point = provide_inlet(vmtk_cap_polydata(polydata))
+        inlet_ind = np.argmin([get_distance(inlet_point, p) for p in center])
+    else:
+        inlet_ind = area.index(max(area))
+
     # Store the center and area
-    inlet_ind = area.index(max(area))
     if case_path is not None:
         info = {"inlet": center[inlet_ind].tolist(), "inlet_area": area[inlet_ind]}
         p = 0
@@ -117,10 +124,34 @@ def compute_centers(polydata, case_path=None):
     return inlet_center, center_
 
 
+def provide_inlet(surface):
+    """
+    Get inlet from user selected point on an input surface
+
+    Args:
+        surface (vtkPolyData): Surface model.
+
+    Returns:
+        points (list): Point located on the inlet, as list
+
+    """
+    print("-- Please select the boundary representing the inlet in the interactive window.")
+    seed_selector = vmtkPickPointSeedSelector()
+    seed_selector.SetSurface(surface)
+    seed_selector.text = "Please select the inlet, \'u\' to undo\n"
+    seed_selector.Execute()
+
+    point_seed_ids = seed_selector.GetTargetSeedIds()
+    get_point = surface.GetPoints().GetPoint
+    inlet_point = list(get_point(point_seed_ids.GetId(0)))
+
+    return inlet_point
+
+
 def provide_relevant_outlets(surface, dir_path=None):
     """
     Get relevant outlets from user
-    selected points on a input surface.
+    selected points on an input surface.
 
     Args:
         surface (vtkPolyData): Surface model.
@@ -154,21 +185,22 @@ def provide_relevant_outlets(surface, dir_path=None):
     return points
 
 
-def get_inlet_and_outlet_centers(surface, base_path, flowext=False):
+def get_inlet_and_outlet_centers(surface, base_path, flowext=False, select_inlet=False):
     """Get the centers of the inlet and outlets.
 
     Args:
         surface (vtkPolyData): An open surface.
         base_path (str): Path to the case file.
         flowext (bool): Turn on/off flow extension.
+        select_inlet (bool): Let user manually select inlet if true
 
     Returns:
         inlet (list): A flatt list with the point of the inlet
         outlet (list): A flatt list with the points of all the outlets.
     """
     # Check if info exists
     if flowext or not path.isfile(base_path + "_info.json"):
-        compute_centers(surface, base_path)
+        compute_centers(surface, base_path, select_inlet)
 
     # Open info
     parameters = get_parameters(base_path)
@@ -187,7 +219,7 @@ def get_inlet_and_outlet_centers(surface, base_path, flowext=False):
             outlets += parameters["outlet%d" % i]
 
     if inlet == [] and outlets == []:
-        inlet, outlets = compute_centers(surface, base_path)
+        inlet, outlets = compute_centers(surface, base_path, select_inlet)
 
     return inlet, outlets
 
@@ -401,7 +433,7 @@ def check_if_surface_is_merged(surface, centerlines, output_filepath):
     for i in range(centerlines.GetNumberOfLines()):
         line_to_compare = vmtk_resample_centerline(lines_to_compare[i], length=0.1)
         line_to_check = vmtk_resample_centerline(extract_single_line(lines_to_check, i), length=0.1)
-        # Compare distance between points along both centerliens
+        # Compare distance between points along both centerlines
         n = min([line_to_check.GetNumberOfPoints(), line_to_compare.GetNumberOfPoints()])
         tolerance = get_centerline_tolerance(line_to_compare) * 500
         for j in range(n):
@@ -417,16 +449,15 @@ def check_if_surface_is_merged(surface, centerlines, output_filepath):
                                     " poly_ball_size.").format(tmp_path))
 
 
-def prepare_output_surface(surface, original_surface, new_centerline, output_filepath,
-                           test_merge=False, changed=False, old_centerline=None,
-                           removed=[[1e9, 1e9, 1e9]]):
-    """After manipulation preparing the surface for output. This method clipps the
-    outlets, slightly smooths the surface, and (potentially) tests if the surface is is
+def prepare_output_surface(surface, original_surface, new_centerline, output_filepath, test_merge=False, changed=False,
+                           old_centerline=None, removed=[[1e9, 1e9, 1e9]]):
+    """After manipulation preparing the surface for output. This method clips the
+    outlets, slightly smooths the surface, and (potentially) tests if the surface is
     merged.
 
     Args:
         surface (vtkPolyData): The new surface after manipulation.
-        original_surface (vtkPolyData): The original surface inputed for manipulation.
+        original_surface (vtkPolyData): The original surface input for manipulation.
         new_centerline (vtkPolyData): The centerline after manipulation.
         output_filepath (str): The user-defined path to the output.
         test_merge (bool): Turn on/off testing if the surface is merged.
@@ -553,7 +584,7 @@ def prepare_output_surface(surface, original_surface, new_centerline, output_fil
 
 
 def attach_clipped_regions_to_surface(surface, clipped, center):
-    """Check the connectivty of a clipped surface, and attach all sections which are not
+    """Check the connectivity of a clipped surface, and attach all sections which are not
     closest to the center of the clipping plane.
 
     Args:
@@ -588,23 +619,21 @@ def attach_clipped_regions_to_surface(surface, clipped, center):
 
 
 def prepare_voronoi_diagram(capped_surface, centerlines, base_path, smooth, smooth_factor, no_smooth, no_smooth_point,
-                            voronoi, pole_ids, resampling_length, absolute=False, upper=None):
+                            voronoi, pole_ids, resampling_length):
     """
     Compute and smooth voronoi diagram of surface model.
 
     Args:
-        capped_surface (polydata): Cappedsurface model to create a Voronoi diagram of.
+        capped_surface (polydata): Capped surface model to create a Voronoi diagram of.
         base_path (str): Absolute path to surface model path.
         voronoi (vtkPolyData): Voronoi diagram.
         pole_ids (vtkIDList): Pole ids of Voronoi diagram.
         smooth (bool): Voronoi is smoothed if True.
         smooth_factor (float): Smoothing factor for voronoi smoothing.
         centerlines (vtkPolyData): Centerlines throughout geometry.
         no_smooth (bool): Part of Voronoi is not smoothed.
-        no_smooth_point (vtkPolyData): Point which defines unsmoothed area.
+        no_smooth_point (vtkPolyData): Point which defines un-smoothed area.
         resampling_length (float): Length of resampling the centerline.
-        absolute (bool): Turn on/off absolute values for the smoothing. Default is off.
-        upper (int): Set an upper limit for the smoothing factor. Default is None.
 
     Returns:
         voronoi (vtkPolyData): Voronoi diagram of surface.
@@ -694,7 +723,7 @@ def compute_centerlines(inlet, outlet, filepath, surface, resampling=1.0, smooth
 def prepare_surface(base_path, surface_path):
     """
     Clean and check connectivity of surface.
-    Capps or uncapps surface at inlet and outlets.
+    Caps or uncaps surface at inlet and outlets.
 
     Args:
         base_path (str): Absolute path to base folder.
@@ -810,13 +839,13 @@ def get_no_smooth_cl(capped_surface, centerlines, base_path, smooth, no_smooth,
     """
     Extract a section where the Voronoi should not be smoothed
      Args:
-        capped_surface (polydata): Cappedsurface model to create a Voronoi diagram of.
+        capped_surface (polydata): Capped surface model to create a Voronoi diagram of.
         centerlines (vtkPolyData): Centerlines throughout geometry.
         base_path (str): Absolute path to surface model path.
         smooth (bool): Voronoi is smoothed if True.
         no_smooth (bool): Part of Voronoi is not smoothed.
         voronoi (vtkPolyData): Voronoi diagram.
-        no_smooth_point (vtkPolyData): Point which defines unsmoothed area.
+        no_smooth_point (vtkPolyData): Point which defines un-smoothed area.
         pole_ids (vtkIDList): Pole ids of Voronoi diagram.
         resampling_length (float): Length of resampling the centerline.
         region_points (list): Flatten list with the region points
@@ -896,7 +925,5 @@ def get_no_smooth_cl(capped_surface, centerlines, base_path, smooth, no_smooth,
 
     no_smooth_cl = vtk_merge_polydata(no_smooth_segments)
     write_polydata(no_smooth_cl, no_smooth_path)
-    # else:
-    #    no_smooth_cl = read_polydata(no_smooth_path)
 
     return no_smooth_cl