add fiber path calculation

models/view_cylinder2.py

@@ -4,6 +4,9 @@
 from boneload import *
 
 class Plane:
+    """construct a plane mesh and create associated mapper & actor
+    (used to see the plane orientation - not used for calculations) 
+    """
     def __init__(self, c, p1, p2):
         source = vtk.vtkPlaneSource()
         source.SetOrigin(c.x[0], c.x[1], c.x[2])
@@ -23,12 +26,42 @@ class Plane:
         self.actor = actor
 
 
+class ClipPolyData:
+    """cut polydata with a plane in 2 sets of parts and keep 1 set
+    input: point and normal vector defining the plane
+    """
+    def __init__(self, polydata, point, normal):
+        plane = vtk.vtkPlane()
+        plane.SetOrigin(point.x[0],point.x[1],point.x[2] )
+        plane.SetNormal(normal.x[0],normal.x[1],normal.x[2])
+
+        clip = vtk.vtkClipPolyData()
+        clip.SetClipFunction(plane)
+        clip.SetInputData(polydata)
+        clip.Update()
+
+        mapper = vtk.vtkPolyDataMapper()
+        mapper.SetInputConnection( clip.GetOutputPort() )
+
+        actor = vtk.vtkActor()
+        actor.SetMapper(mapper)        
+        actor.GetProperty().SetColor( colors.GetColor3d('orange') )
+
+        self.plane = plane
+        self.clip = clip
+        self.mapper = mapper
+        self.actor = actor
+
+
 class PlaneCut:
-    def __init__(self, polydata, c, n):
+    """cut a polydata with a plane and keep the intersection line
+    input: point and normal vector defining the plane
+    """
+    def __init__(self, polydata, point, normal):
 
         plane = vtk.vtkPlane()
-        plane.SetOrigin(c.x[0], c.x[1], c.x[2])
-        plane.SetNormal(n.x[0], n.x[1], n.x[2])
+        plane.SetOrigin(point.x[0], point.x[1], point.x[2])
+        plane.SetNormal(normal.x[0], normal.x[1], normal.x[2])
 
         cutter = vtk.vtkCutter()
         cutter.SetCutFunction(plane)
@@ -50,11 +83,14 @@ class PlaneCut:
 
 
 class ClosestPart:
+    """keep the part of a polydata which is closest to a given point
+    """
     def __init__(self, outputport, point):
         cfilter = vtk.vtkPolyDataConnectivityFilter() 
         cfilter.SetExtractionModeToClosestPointRegion()	
         cfilter.SetClosestPoint(point.x[0], point.x[1], point.x[2]) 
         cfilter.SetInputConnection(outputport)
+        cfilter.Update()
 
         mapper = vtk.vtkPolyDataMapper()
         mapper.SetInputConnection( cfilter.GetOutputPort())
@@ -69,6 +105,42 @@ class ClosestPart:
         self.actor = actor
 
 
+def get_segments(poly):
+    """extract points coordinates and connectivities of the segments of
+    a polydata made of a series of line segments
+    """
+    lines = poly.GetLines()
+    points = poly.GetPoints()
+    # print(f'nb of segments= {lines.GetNumberOfCells()}')
+
+    # fill points
+    pts = []
+    for i in range(points.GetNumberOfPoints()):
+        pts.append(Pt(points.GetPoint(i)))
+    # print(f'pts={[ str(p) for p in pts]}')
+
+    ids = []
+    idList = vtk.vtkIdList()
+    lines.InitTraversal()
+    while lines.GetNextCell(idList):
+        npts = idList.GetNumberOfIds()
+        if npts!=2:
+            print("intersection contains bad cells")
+            continue
+        ids.append(( idList.GetId(0),idList.GetId(1) ))
+    return pts, ids
+
+
+def compute_length(pts, ids):
+    """compute the total length of a series of line segments
+    """
+    length=0.
+    for v in ids:
+        length+=abs(pts[v[1]]-pts[v[0]])
+    return length
+
+
+
 if __name__=="__main__":
 
     # load meshes
@@ -101,35 +173,76 @@ if __name__=="__main__":
         
         force = (area/3)*direct      # compute nodal force
 
-        plane = Plane(centre, centre+nunit, targetP) 
 
+        plane = Plane(centre, centre+nunit, targetP) 
         nclip = -(targetP - (centre+(targetP*nunit)*nunit)).normalized()
-
-        # clip = vtk.vtkPolyDataPlaneClipper() # vtk 9..
-
-        cplane = vtk.vtkPlane()
-        cplane.SetOrigin(centre.x[0],centre.x[1],centre.x[2] )
-        cplane.SetNormal(nclip.x[0],nclip.x[1],nclip.x[2])
-        clip = vtk.vtkClipPolyData()
-        clip.SetClipFunction(cplane)
-        clip.SetInputData(muscle.polydata)
-        clip.Update()
-
-        mapper = vtk.vtkPolyDataMapper()
-        mapper.SetInputConnection( clip.GetOutputPort())
-
-        actor = vtk.vtkActor()
-        actor.SetMapper(mapper)        
-        actor.GetProperty().SetColor( colors.GetColor3d('orange') )
-        # actor.GetProperty().SetLineWidth(2)
-
-
-
-
-        planecut = PlaneCut(clip.GetOutput(), centre, nunit.cross(targetP-centre))
+        planepart = ClipPolyData(muscle.polydata, centre, nclip)
+        planecut = PlaneCut(planepart.clip.GetOutput(), centre, nunit.cross(targetP-centre))
         closest = ClosestPart(planecut.cutter.GetOutputPort(), centre)
         
+        # compute "s(r)"
+        pts, ids = get_segments(closest.cfilter.GetOutput())
+        s = compute_length(pts, ids)
+        print(f's(r)={s}')
+
+        ed_by_nods = [ [] for p in pts ]
+        # sort segments
+        for v in ids:
+            ed_by_nods[ v[0] ].append( v )
+            ed_by_nods[ v[1] ].append( v )
+        print(f'ed_by_nods={ed_by_nods}')    
+
+        # look for first segment
+        for curp,eds in enumerate(ed_by_nods):
+            if len(eds)==1:
+                edge = eds[0]
+                nextp = edge[0]
+                if nextp==curp:
+                    nextp=edge[1]
+                break
+
+        edges = []
+        if edge[0] == curp:
+            edges.append( edge )
+        else:
+            edges.append( (edge[1], edge[0]) )
+
+        # print(f'edge={edge}')
+        # print(f'nextp={nextp}')
+
+        while True:
+            curp = nextp
+            eds = ed_by_nods[curp]
+            # print(f'eds={eds}')
+            if len(eds)==1:
+                # print('last one!')
+                edge = eds[0]
+                if edge[0] == curp:
+                    edges.append( edge )
+                else:
+                    edges.append( (edge[1], edge[0]) )
+                break
+            lastedge = edge
+            edge = eds[0]
+            if edge == lastedge:
+                edge = eds[1]
+            
+            if edge[0] == curp:
+                edges.append( edge )
+            else:
+                edges.append( (edge[1], edge[0]) )
+
+            nextp = edge[0]
+            if nextp==curp:
+                nextp=edge[1]
+
+            # print(f'edge={edge}')
+            # print(f'nextp={nextp}') 
+            # input()      
+
+        print(f'edges={edges}')
 
+    
 
     # -- view
 
@@ -139,8 +252,7 @@ if __name__=="__main__":
     view.addActors(muscle.actor)  
     view.addActors(muscleF.actor)
     view.addActors(plane.actor)  
-    view.addActors(actor) 
-    
+    view.addActors(planepart.actor) 
     view.addActors(planecut.actor) 
     view.addActors(closest.actor)