Skip to content
Snippets Groups Projects
Commit da298208 authored by Boman Romain's avatar Boman Romain
Browse files

accurate curvature

parent 9a9ae522
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
models/view_cylinder2.py
+ 128
− 102
View file @ da298208
......@@ -140,7 +140,7 @@ def compute_length(pts, ids):
for v in ids:
dl = abs(pts[v[1]]-pts[v[0]])
length+=dl
print(f'\tdl={dl}')
# print(f'\tdl={dl}')
return length
......@@ -221,6 +221,101 @@ def sort_segments(pts, ids):
return sorted_ids
def sort_vertices(pts, ids, centre):
"""build a list of vertex indices of an open piecewise linear path.
the points are sorted by increasing distance to "centre" (Pt object)
input: pts: a list of points coordinates (Pt objects)
ids: a list of tuples of indexing the 2 vertices of each segment
"""
# print(f'ids={ids}')
sorted_ids = sort_segments(pts, ids)
# print(f'sorted_ids={sorted_ids}')
# check whether "centre" is the last vertex
# reverse the list if it is the case
if( abs(centre-pts[sorted_ids[-1][0]]) < abs(centre-pts[sorted_ids[0][0]]) ):
sorted_ids = [ (seg[1], seg[0]) for seg in reversed(sorted_ids)]
# print(f'sorted_ids={sorted_ids}')
# build sorted vertex indices
sorted_v = [ seg[0] for seg in sorted_ids ] # append first vertex of each segment
sorted_v.append( sorted_ids[-1][1] ) # append last vertex
# print(f'sorted_v={sorted_v}')
return sorted_v
def build_clean_path(pts, sorted_v):
# build a sorted list of coordinates along the muscle fibre
fpath = []
for iv in sorted_v:
fpath.append( pts[ iv ] )
# longest segment
dlmax = 0.0
for i in range(len(fpath)-1):
dl = abs(fpath[i+1]-fpath[i])
dlmax = max(dlmax, dl)
# print(f'dlmax={dlmax}')
# remove colinear segments or too small segments
while True:
to_remove=None
for i in range(1, len(fpath)-1):
p0 = fpath[i-1]
p1 = fpath[i]
p2 = fpath[i+1]
d1 = (p1-p0)
d2 = (p2-p1)
d1u = d1.normalized()
d2u = d2.normalized()
# print(f'd1u.cross(d2u)={abs(d1u.cross(d2u))}')
if(abs(d1u.cross(d2u))<1.e-2):
to_remove = i
break
if abs(d1)<1.e-3*dlmax or abs(d1)<1.e-3*dlmax:
to_remove = i
break
if to_remove!=None:
# print(f'removing vertex at {fpath[to_remove]}')
del fpath[to_remove]
else:
break
return fpath
def build_local_axes(fpath, nunit):
p0 = fpath[0] # origin
p1 = fpath[1] # next point
xaxis = (p1-p0).normalized()
yaxis = nunit
return p0, xaxis, yaxis
def compute_curvature(fpath, fpath2, nunit):
# fit a polynomial
fitpts = fpath2[2:0:-1] + fpath[1:3]
# print(f'len(fitpts)={len(fitpts)}')
order = len(fitpts)-1
# build local axes
p0, xaxis, yaxis = build_local_axes(fpath, nunit)
import numpy as np
x = [ (p-p0)*xaxis for p in fitpts ]
y = [ (p-p0)*yaxis for p in fitpts ]
z = np.polyfit(x, y, order)
polynom = np.poly1d(z)
polyx = np.linspace(x[0],x[-1],100) # only for display
polyy = polynom(polyx) # only for display
# https://www.math24.net/curvature-radius
polynomd = polynom.deriv()
polynomdd = polynomd.deriv()
curv = abs(polynomdd(0.0))/math.pow(1.0+polynomd(0.0)**2, 3/2)
return curv, polyx, polyy
if __name__=="__main__":
......@@ -240,8 +335,10 @@ if __name__=="__main__":
targetP = Pt(focalnodes[0] )
total_area=0.
notri = 140
# notri = 141
# notri = 140
# notri = 1
notri = 2
# notri = 3
for tri in [ tris[notri] ]:
p1, p2, p3 = [ Pt(nodes[tri[i]]) for i in range(3) ]
......@@ -264,118 +361,44 @@ if __name__=="__main__":
planepart = ClipPolyData(muscle.polydata, centre, nclip)
planecut = PlaneCut(planepart.clip.GetOutput(), centre, nunit.cross(targetP-centre))
closest = ClosestPart(planecut.cutter.GetOutputPort(), centre)
pts, ids = get_segments(closest.cfilter.GetOutput())
# compute "s(r)"
pts, ids = get_segments(closest.cfilter.GetOutput())
s = compute_length(pts, ids)
s = compute_length(pts, ids)
print(f's(r)={s}')
print(f'ids={ids}')
sorted_ids = sort_segments(pts, ids)
print(f'sorted_ids={sorted_ids}')
# check whether "centre" is the last vertex
# reverse the list if it is the case
if( abs(centre-pts[sorted_ids[-1][0]]) < 1.0e-8 ):
sorted_ids = [ (seg[1], seg[0]) for seg in reversed(sorted_ids)]
print(f'sorted_ids={sorted_ids}')
# build sorted vertex indices
sorted_v = [ seg[0] for seg in sorted_ids ] # append first vertex of each segment
sorted_v.append( sorted_ids[-1][1] ) # append last vertex
print(f'sorted_v={sorted_v}')
# | t_1-t_0 |
# curvature = lim -------------
# p_1->p_0 | p_1-p_0 |
# t = tangent evaluated at point p
# tangents = []
# centres = []
# for i,seg in enumerate(sorted_ids):
# p0 = pts[ seg[0] ]
# p1 = pts[ seg[1] ]
# tangents.append( (p1-p0).normalized() )
# # first seg: use the first vertex as centre because the face was cut
# if i==0:
# c=p0
# else:
# c=(p1+p0)/2
# centres.append(c)
# curv=0.
# if (len(centres)>1):
# t0= tangents[0]
# t1= tangents[1]
# p0= centres[0]
# p1= centres[1]
# curv = abs(t1-t0)/abs(p1-p0)
# print(f'curv={curv}')
# build local axes
seg = sorted_ids[0]
p0 = pts[ seg[0] ]
p1 = pts[ seg[1] ]
xaxis = (p1-p0).normalized()
yaxis = n
# build a sorted list of coordinates along the msucle fibre
fpath = []
for iv in sorted_v:
fpath.append( pts[ iv ] )
# longest segment
dlmax = 0.0
for i in range(len(fpath)-1):
dl = abs(fpath[i+1]-fpath[i])
dlmax = max(dlmax, dl)
print(f'dlmax={dlmax}')
# remove colinear segments or too small segments
while True:
to_remove=None
for i in range(1, len(fpath)-1):
p0 = fpath[i-1]
p1 = fpath[i]
p2 = fpath[i+1]
d1 = (p1-p0)
d2 = (p2-p1)
d1u = d1.normalized()
d2u = d2.normalized()
# print(f'd1u.cross(d2u)={abs(d1u.cross(d2u))}')
if(abs(d1u.cross(d2u))<1.e-3):
to_remove = i
break
if abs(d1)<1.e-3*dlmax or abs(d1)<1.e-3*dlmax:
to_remove = i
break
if to_remove!=None:
print(f'removing vertex at {fpath[to_remove]}')
del fpath[to_remove]
else:
break
# append symetrical point before the first one
# p0 = pts[ seg[0] ]
# p1 = pts[ seg[1] ]
# fpath = [ p0 - (p1-p0) ] + fpath
# project vertices onto the plane (display)
xpath = []
ypath = []
for p in fpath:
xpath.append((p-p0)*xaxis)
ypath.append((p-p0)*yaxis)
# front part of the cut
planepart2 = ClipPolyData(muscle.polydata, centre, -nclip)
planecut2 = PlaneCut(planepart2.clip.GetOutput(), centre, nunit.cross(targetP-centre))
closest2 = ClosestPart(planecut2.cutter.GetOutputPort(), centre)
pts2, ids2 = get_segments(closest2.cfilter.GetOutput())
sorted_v = sort_vertices(pts, ids, centre)
sorted_v2 = sort_vertices(pts2, ids2, centre)
fpath = build_clean_path(pts, sorted_v)
fpath2 = build_clean_path(pts2, sorted_v2)
curv, polyx, polyy = compute_curvature(fpath, fpath2, nunit)
print(f'curv={curv}')
# display things
# project vertices of fpath onto the plane for display
p0, xaxis, yaxis = build_local_axes(fpath, nunit)
xpath = [ (p-p0)*xaxis for p in reversed(fpath2) ] + [ (p-p0)*xaxis for p in fpath ]
ypath = [ (p-p0)*yaxis for p in reversed(fpath2) ] + [ (p-p0)*yaxis for p in fpath ]
import matplotlib.pyplot as plt
import numpy as np
plt.plot(xpath, ypath, 'o-', label='path')
plt.plot(polyx, polyy, 'r-', label='polynomial')
plt.xlabel('x')
plt.ylabel('y')
plt.title('muscle fiber shape')
plt.grid(True)
plt.legend()
#plt.savefig("test.png")
plt.show()
# -- view
......@@ -389,6 +412,7 @@ if __name__=="__main__":
view.addActors(planepart.actor)
view.addActors(planecut.actor)
view.addActors(closest.actor)
view.addActors(closest2.actor)
surface.actor.GetProperty().SetOpacity(0.1)
......@@ -402,5 +426,7 @@ if __name__=="__main__":
muscleF.actor.GetProperty().SetPointSize(15.0)
muscleF.actor.GetProperty().SetColor( colors.GetColor3d('green') )
closest2.actor.GetProperty().SetColor( colors.GetColor3d('pink') )
view.interact()
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment