"""
Practical functionality to process rsml-formated mtg
"""
[docs]
def plant_vertices(g):
""" return the list of mtg vertices that represent plants """
return g.vertices(scale=1)
[docs]
def root_vertices(g):
""" return the list of mtg vertices that represent root axes """
return g.vertices(scale=g.max_scale())
[docs]
def root_tree(g, suborder=None):
""" return the list of root axes in topological order
If suborder is given, it should be a dictionary of (root-id,value) which is
used to sort sibling root w.r.t. their respective value.
"""
if suborder is None:
sort = lambda x:x[::-1]
else:
sort = lambda x: sorted(x, key=suborder.get)[::-1]
# init with axes with no parent
axes = sort([a for a in root_vertices(g) if g.parent(a) is None])
# parse root axes in depth-first-order
tree = []
while len(axes):
axe = axes.pop()
tree.append(axe)
axes.extend(sort(g.children(axe)))
return tree
[docs]
def root_order(g, tree=None):
""" return a dictionary of the (numeric) axe order
The order is select as:
- the value of the 'order' property, if present
- otherwise, 1 for axe with no parent or the parent order +1
tree is the optional list of root id in `g`. If not given, it is computed
"""
# init order dict
if 'order' in g.property_names():
order = g.property('order').copy()
else:
order = {}
if tree is None:
tree = root_tree(g)
# parse all axes s.t. parent are processed before children
for root in tree:
parent = g.parent(root)
order.setdefault(root, 1 if parent is None else order[parent]+1)
return order
[docs]
def hausdorff_distance(polyline1,polyline2):
"""
Compute the hausdorff distance from `polyline1` to `polyline2`
:Inputs:
`polyline1`:
a (k,n1) array for the n1 points of the 1st polyline in k-dimension
`polyline2`:
a (k,n2) array for the n2 points of the 2nd polyline in k-dimension
:Output:
The hausdorff distance:
max( D(polyline1,polyline2), D(polyline2,polyline1) )
where
D(p1,p2) = max_(i in p1) |p1[i] - closest-projection-on-p2|
"""
import numpy as _np
p1 = _np.asfarray(polyline1)
p2 = _np.asfarray(polyline2)
norm = lambda x: (x**2).sum(axis=0)**.5
def max_min_dist(points, polyline):
v1 = polyline[:,:-1] # 1st segment vertex, shape (k,n2-1)
v2 = polyline[:, 1:] # 2nd segment vertex, shape (k,n2-1)
sdir = v2-v1 # direction vector of segment
lsl = norm(sdir) # distance between v1 and v2
lsl = _np.maximum(lsl,2**-5)
sdir /= lsl # make sdir unit vectors
# distance from v1 to the projection of points on segments
# disallow projection out of segment: values are in [0,lsl]
on_edge = ((points[:,:,None]-v1[:,None,:])*sdir[:,None,:]).sum(axis=0) # (n1,n2-1)
on_edge = _np.minimum(_np.maximum(on_edge,0),lsl[None,:])
# points projection on sdir
nproj = v1[:,None,:] + on_edge[None,:,:]*sdir[:,None,:] # (k,n1,n2-1)
# distance from points to "points projection on sdir"
return norm(nproj - points[:,:,None]).min(axis=1).max()
return max(max_min_dist(p1,p2), max_min_dist(p2,p1))
