# -*- python -*-
#
# OpenAlea.mtg
#
# Copyright 2008-2012 INRIA - CIRAD - INRA
#
# File author(s): Christophe Pradal <christophe.pradal.at.cirad.fr>
#
# Distributed under the Cecill-C License.
# See accompanying file LICENSE.txt or copy at
# http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html
#
# OpenAlea WebSite : http://openalea.gforge.inria.fr
#
################################################################################
''' This module implement a solver to build 3D representation of a plant mockup
based on various infomation.
:Principles:
Use topologic information and associated properties to reconstruct a 3D representation
of a plant. The 3D representation must satisfy all the intra and inter topological constraints
and the properties defined by the user.
:Algorithm:
:Examples:
.. todo::
- define a set of algorithms that can be used independently (as plug'in).
- for each algo:
* Decompose the problem into subproblems by computing the frontere.
=>define a color based on the properties and the semantic.
=> define a partial order for solving each problems
* Solve each subproblems and assemble the solutions
* Return the solution as properties
- define a generic method (or methods) to build a geometric object for computed values and properties.
- diameter
- length vs points
- angles
- remove noise from points
- translate the branch at the extremities of the cylinder.
'''
from math import sqrt, pi
from .. import traversal
from .dresser import DressingData
from .. import algo
from openalea.mtg.mtg import colored_tree, PropertyTree
from openalea.mtg.plantframe import turtle
from openalea.plantgl.all import *
error = True
epsilon = 1e-5
[docs]
class PlantFrame(object):
''' Engine to compute the geometry of a plant based on
its topological description and parameters.
'''
def __init__(self, g, *args, **kwds):
''' Compute a geometric representation of a plant.
Compute the skeleton of the MTG with a Frame linked to each elements.
:Parameters:
- g: MTG representing the plant architecture at different scales.
'''
# Express everything in degrees or radians.
# Do not mix both.
# Define the inference algorithms to compute missing parameters.
# Add new algorithms as:
# - Visitors
# - plug'in.
# - derivation
# 1. Get the different parameters
# 2. Build for each element at a given scale the frames.
# 3 Solve the different constraints by using a plug'in mechanism
# 4. Design new algo and define new parameters.
# 5. Use WeberPenn with MTG
self.g = g
self.max_scale = g.max_scale()
# Global parameters
self.origin = kwds.get('origin', (0,0,0))
self.dresser = kwds.get('DressingData', DressingData())
# properties defined for each vertex.
# length define the length of a vertex
self.length = self._extract_properties('Length', kwds)
self.top_diameter = self._extract_properties('TopDiameter', kwds)
self.bottom_diameter = self._extract_properties('BottomDiameter', kwds)
self.diameter = self.top_diameter
# Absolute coordinate of the components
self.xx= self._extract_properties('XX', kwds)
self.yy = self._extract_properties('YY', kwds)
self.zz = self._extract_properties('ZZ', kwds)
# Absolute euler angles in degree
self.aa = self._extract_properties('AA', kwds)
self.bb = self._extract_properties('BB', kwds)
self.cc = self._extract_properties('CC', kwds)
self.alpha = self._extract_properties('Alpha', kwds)
# Relative angles: TO BE DEFINED: alpha, beta, rolll
# Curves associated to an element
# Geometry associated to an element with mtg parameters
# User define axes: Return True or False
self.new_axe = self._extract_properties('Axe', kwds)
if not self.new_axe:
self.new_axe = dict((vid, True) for vid in g.vertices_iter() if g.edge_type(vid) == '+')
self.axes = {}
# Exclude some elements depending on their class
self.exclude = kwds.get('Exclude', [])
self._length = {}
self._diameter = {}
self._surface = {}
self._volume = {}
self._segmentvec = {}
self._inertia = {}
self._compute_global_data()
self.propagate_constraints()
def _extract_properties(self, name, kwds):
''' Extract the property from properties of the mtg
or from a user given function.
'''
d = {}
func = kwds.get(name)
if callable(func):
# Compute the value for all vertices
all_values = ((vid,func(vid)) for vid in self.g.vertices_iter())
# Select only those which are defined
values = ( (vid, fvid) for vid, fvid in all_values if fvid is not None)
d = dict(values)
elif isinstance(func, str):
d = self.g.property(func)
name_property = self.g.property(name).copy()
name_property.update(d)
if name in ['XX', 'YY', 'ZZ', 'Length']:
factor = 1./self.dresser.length_unit
elif name in ['TopDiameter', 'BottomDiameter']:
factor = 1./self.dresser.diameter_unit
elif name in ['Alpha', 'AA', 'BB', 'CC']:
factor = 1./self.dresser.alpha_unit
else:
return name_property
for k, v in name_property.items():
name_property[k] = factor * v
return name_property
def _compute_global_data(self):
"""
Scale the values according to units.
"""
g = self.g
self.points = {}
points = self.points
for vid in self.xx:
try:
points[vid] = self.xx[vid], self.yy[vid], self.zz[vid]
except:
continue
self.euler_angles = {}
angles = self.euler_angles
for vid in self.aa:
try:
angles[vid] = self.aa[vid], self.bb[vid], self.cc[vid]
except:
continue
if not self.new_axe:
self.new_axe = [v for v in g.vertices_iter(scale = self.max_scale) if (g.edge_type(v) == '+' or g.parent(v) is None) or g.class_name(v) in self.exclude]
# Method to compute the axes and their order.
self._compute_axes()
def _compute_axes(self):
"""
Compute the axes with their order.
"""
if self.axes:
return
g = self.g
max_scale = self.max_scale
marked = {}
axes = self.axes
new_axe = self.new_axe
exclude = self.exclude
self.max_order_axes = 0
def ancestors(v):
while v is not None:
yield v
marked[v] = True
if v in new_axe:
break
v = g.parent(v)
def order(v):
_order = 0
while v is not None:
if v in new_axe:
_order += 1
v = g.parent(v)
return _order
for root in g.roots_iter(scale=max_scale):
for vid in traversal.post_order(g, root):
if vid in marked or g.class_name(vid) in exclude:
continue
_axe = list(reversed(list(ancestors(vid))))
_order = order(_axe[0])
axes.setdefault(_order, []).append(_axe)
def _get_origin(self, vid):
''' Compute the origin for the vertex `vid`.
'''
# 1. compute the origin of the tree
# Check if a complex has origin
origin = self.origin
vid = self.g.complex(vid)
while vid:
if vid in self.points:
origin = self.points[vid]
break
else:
vid = self.g.complex(vid)
self.origin = origin
return self.origin
[docs]
def run(self, scale= -1):
'''Compute the geometry of the plant.
'''
if scale == -1:
scale = self.g.max_scale()
# TODO: Compute the plantframe for several plants.
root = self.g.roots(scale=scale)[0]
# 1. compute the origin of the tree
# Check if a complex has origin
self.origin = self._get_origin(root)
axes = compute_axes(self.g, root, self.points, self.origin)
# Compute the points
# 1. compute fixed points
#
# Compute diameters
self.algo_diameter()
[docs]
def propagate_constraints(self):
""" Propagate the properties into the whole MTG.
"""
g = self.g
# diameter
max_scale = g.max_scale()
# Between scale constraints
# Copy of the keys because we modify the
for vid in list(self.bottom_diameter.keys()):
scale = g.scale(vid)
if scale == max_scale:
continue
for s in range(scale+1, max_scale+1):
roots = g.component_roots_at_scale_iter(vid, s)
try:
component_id = next(roots)
self.bottom_diameter[component_id] = self.bottom_diameter[vid]
except StopIteration:
pass
# Within scale constraint
# BottomDiameter(x) = TopDiameter(Parent(x))
# WARNING : This is NOT TRUE:
# - You may have several bot_dia values for the same parent.
edge_type = g.property('edge_type')
for vid in self.bottom_diameter:
pid = g.parent(vid)
if pid is not None :
if pid not in self.top_diameter \
and (g.nb_children(pid) == 1 or edge_type.get(vid) == '<'):
self.top_diameter[pid] = self.bottom_diameter[vid]
# WARNING: We may have several extremities for the components.
for vid in list(self.top_diameter.keys()):
scale = g.scale(vid)
if scale == max_scale:
continue
last = None
if self.is_linear(g, vid):
last = self._last_component(g, vid)
if last not in self.top_diameter:
self.top_diameter[last] = self.top_diameter[vid]
for s in range(scale+2, max_scale):
if self.is_linear(g, last):
last = self._last_component(g, last)
if last not in self.top_diameter:
self.top_diameter[last] = self.top_diameter[vid]
################################################################
# Methods that extend MTG
@staticmethod
def _first_component(g, vid):
return next(self.g.component_roots_iter(vid))
@staticmethod
def _last_component(g, vid):
leaves = algo.extremities(g, vid, scale=g.scale(vid)+1, ContainedIn=vid)
return next(leaves)
@staticmethod
def _first_component_at_scale(g, vid, scale):
return next(g.component_roots_at_scale_iter(vid, scale))
@staticmethod
def _last_component_at_scale(g, vid, scale):
leaves = algo.extremities(g, vid, scale=scale, ContainedIn=vid)
return next(leaves)
[docs]
@staticmethod
def is_linear(g, cid):
"""
A complex is linear iff there is only there is only one extremity in it.
"""
leaves = algo.extremities(g, cid, scale=g.scale(cid)+1, ContainedIn=cid)
return len(list(leaves)) == 1
[docs]
@staticmethod
def strahler_order(g, vid):
strahler = {}
for v in traversal.post_order(g, vid):
children_strahler = [strahler[c] for c in g.children(v)]
if children_strahler:
m, M = min(children_strahler), max(children_strahler)
strahler[v] = M if m != M else M+1
else:
strahler[v] = 1
return strahler
###########################################################################################
# Algorithms for the different versions
###########################################################################################
################################################################
# Write the different plugins
################################################################
# Diameter algorithms
[docs]
def algo_diameter(self, mode=1, power = 2):
"""
Compute the radius for each vertices.
Cases:
1. No radius values at all : pipe model
2. Linear interpolation of radius on axes.
3. Pipe model
"""
if not self.top_diameter and not self.bottom_diameter:
mode = 0
if mode == 0:
diameter = self.default_algo_diameter(power)
else:
#diameter = self.advanced_algo_diameter(power)
diameter = self.advanced_algo_diameter2(power)
self._diameter = diameter
return diameter
[docs]
def default_algo_diameter(self, power):
"""
Compute the pipe model on the entire mtg.
When a node has no values, take the deafult value.
"""
g = self.g
max_scale = g.max_scale()
v = next(g.roots_iter(scale=max_scale))
# Compute default diameter
dresser = self.dresser
default_diameter = 1 if not dresser.min_topdia else min(dresser.min_topdia.values())
default_diameter *= 1./self.dresser.diameter_unit
strands = {}
for vid in traversal.post_order(g, v):
strands[vid] = max(sum([strands[c] for c in g.children(vid)]), 1)
diameters = {}
for vid, s in strands.items():
diameters[vid] = default_diameter * pow(s, 1./power)
return diameters
[docs]
def advanced_algo_diameter(self, power, default_diameter=None):
"""
Compute the pipe model on the entire mtg.
There are 4 cases:
for each component:
- traverse in a post_order way all the vertices
- compute strands and diameter
- compute the difference at the root
- compute a diameter value for the strands
- fix problem with wrong power value:
The diameter has to be decreasing.
- if no value for the root, compute the total strands for the tree.
- Then divide by the defined radius values obtain from step 1.
May be improved by a filtering pass.
Select all the segment without ramification (linear) and interpolate the radius.
"""
g, new_map = self.build_mtg_from_radius()
max_scale = g.max_scale()
dresser = self.dresser
if default_diameter is None:
default_diameter = 1 if not dresser.min_topdia else min(dresser.min_topdia.values())
default_diameter *= 1./self.dresser.diameter_unit
default_diameter = default_diameter**power
strands = {}
diameters = {}
# update the defined properties with the new indices
top_diameter = {}
bottom_diameter = {}
for v in g.vertices_iter(scale=2):
old_v = new_map[v]
if old_v in self.top_diameter:
top_diameter[v] = self.top_diameter[old_v]
if old_v in self.bottom_diameter:
bottom_diameter[v] = self.bottom_diameter[old_v]
for k,td in top_diameter.items():
diameters[k] = td**power
###
error_vertex = []
# For each independant sub_systems
for cid in g.vertices_iter(scale=1):
# traverse the tree in a post_order way only for all vid in cid
root = next(g.component_roots_iter(cid))
has_root_diameter = root in diameters or root in bottom_diameter or g.parent(root) in diameters
if has_root_diameter:
pid = g.parent(root)
if pid in diameters:
root_diameter = diameters[g.parent(root)]
else:
root_diameter = max(diameters.get(root, 0),
bottom_diameter.get(root,0)**power)
sorted_vertices = list(traversal.post_order(g, root, complex=cid))
for vid in sorted_vertices:
if vid in diameters:
continue
children = [v for v in g.children(vid) if g.complex(v) == cid]
children_diam = [diameters[v] for v in children if v in diameters]
diam = sum(children_diam)
if has_root_diameter and diam > root_diameter:
if error:
if max(children_diam) < root_diameter:
print("ONE children has a greater radius value than its root.")
print(children, children_diam)
else:
print('WARNING: The pipe model compute at %d for power=%f a too large diameter.'%(vid, power))
print(' -> decrease the power of the pipe model.')
print('root ', root, 'root_diam ', root_diameter, 'current ', diam)
error_vertex.append(new_map[vid])
if diam > 0:
diameters[vid] = diam
strand = sum(strands.get(v,0) for v in children)
if strand > 0:
strands[vid] = strand
elif diam == 0:
strands[vid] = 1
# Solve the boundary condition
# if there are a bottom diam on root
# check
for vid in sorted_vertices:
assert vid in strands or vid in diameters
if has_root_diameter:
#if len(sorted_vertices) == 1:
# diameters[root] = root_diameter
# continue
delta_diameter = root_diameter - diameters.get(root,0)
if delta_diameter < epsilon:
# compute strands default radius
# Add diameter only at the branching where
# Select only the vertices which do not have a diam value.
strands_diameter = default_diameter
vtxs = [v for v in sorted_vertices if v in strands and v not in diameters]
for vid in vtxs:
s = strands[vid]
diameters[vid] = strands_diameter * s
else:
# d**2 = sum(d**2) + sum(strand_diameter**2) = n * cst_strand_diameter**2
n = strands.get(root,0)
if n > 0:
strands_diameter = delta_diameter / n
vtxs = [v for v in sorted_vertices if v in strands]
for vid in vtxs:
diameters[vid] = diameters.get(vid, 0) + strands[vid] * strands_diameter
else:
# Compute the total number of strands for the mtg
nb_leaves = len([v for v in g.vertices_iter(scale=2) if g.is_leaf(v)])
if root not in diameters:
strands_diameter = default_diameter
elif root not in strands:
continue
else:
strands_diameter = diameters[root] / (nb_leaves - strands[root])
vtxs = [v for v in sorted_vertices if v in strands]
for vid in vtxs:
diameters[vid] = diameters.get(vid, 0) + strands[vid] * strands_diameter
# compute the final diameters
factor = 1./power
result = dict(( (new_map[v], d**factor) for v, d in diameters.items()))
self.error_vertex = error_vertex
if error_vertex:
print('Warnings on diameter for %d vertices'%len(error_vertex))
return result
[docs]
def advanced_algo_diameter2(self, power, default_diameter=None):
"""
Compute the pipe model on the entire mtg.
There are 4 cases:
for each component:
- traverse in a post_order way all the vertices
- compute strands and diameter
- compute the difference at the root
- compute a diameter value for the strands
- fix problem with wrong power value:
The diameter has to be decreasing.
- if no value for the root, compute the total strands for the tree.
- Then divide by the defined radius values obtain from step 1.
May be improved by a filtering pass.
Select all the segment without ramification (linear) and interpolate the radius.
"""
trees = self.decompose_radius()
g = self.g
max_scale = g.max_scale()
dresser = self.dresser
exclude = self.exclude
n = 0
if default_diameter is None:
default_diameter = 1 if not dresser.min_topdia else min(dresser.min_topdia.values())
default_diameter *= 1./self.dresser.diameter_unit
default_diameter = default_diameter**power
strands = {}
diameters = {}
# update the defined properties with the new indices
top_diameter = self.top_diameter
bottom_diameter = self.bottom_diameter
for k,td in top_diameter.items():
if g.scale(k) == max_scale:
diameters[k] = td**power
###
error_vertex = []
# For each independant sub_systems
for tree in trees.values():
# traverse the tree in a post_order way only for all vid in cid
root = tree.root
pid = g.parent(root)
has_root_diameter = root in bottom_diameter or pid in diameters or g.parent(pid) in diameters
strands.clear()
if has_root_diameter:
if pid in diameters:
root_diameter = diameters[pid]
else:
root_diameter = max(diameters.get(g.parent(pid), 0),
bottom_diameter.get(root,0)**power)
# Exclude some nodes
for vid in traversal.post_order(tree, root):
assert vid not in diameters
if g.class_name(vid) in exclude:
continue
children_diam = [diameters[v] for v in g.children(vid) if v in diameters and g.class_name(v) not in exclude]
diam = sum(children_diam)
if has_root_diameter and diam > root_diameter:
if error:
if max(children_diam) < root_diameter:
print("ONE children has a greater radius value than its root.")
print(g.children(vid), children_diam)
else:
print('WARNING: The pipe model compute at %d for power=%f a too large diameter.'%(vid, power))
print(' -> decrease the power of the pipe model.')
print('root ', root, 'root_diam ', root_diameter, 'current ', diam)
error_vertex.append(vid)
if diam > 0:
diameters[vid] = diam
strand = sum(strands.get(v,0) for v in g.children(vid) if g.class_name(v) not in exclude)
if strand > 0:
strands[vid] = strand
elif diam == 0:
strands[vid] = 1
# Solve the boundary condition
# if there are a bottom diam on root
if has_root_diameter:
delta_diameter = root_diameter - diameters.get(root,0)
if delta_diameter < epsilon:
# compute strands default radius
# Add diameter only at the branching where
# Select only the vertices which do not have a diam value.
strands_diameter = default_diameter
for vid, nb_strands in strands.items():
if vid not in diameters:
diameters[vid]= nb_strands * strands_diameter
else:
# d**2 = sum(d**2) + sum(strand_diameter**2) = n * cst_strand_diameter**2
n = strands.get(root,0)
if n > 0:
strands_diameter = delta_diameter / n
for vid, nb_strands in strands.items():
diameters[vid] = diameters.get(vid,0) + nb_strands * strands_diameter
else:
n+=1
# Compute the total number of strands for the mtg
if strands.get(root,0)==0:
assert not bool(strands)
assert root in diameters
continue
def special_leaf(v):
return not any(cid for cid in g.children(v) if g.class_name(cid) not in exclude)
if root not in diameters:
strands_diameter = default_diameter
else:
nb_leaves = len([v for v in traversal.pre_order2(g, root) if special_leaf(v)])
strands_diameter = diameters[root] / (nb_leaves - strands[root])
for vid, nb_strands in strands.items():
diameters[vid] = diameters.get(vid,0) + nb_strands * strands_diameter
# compute the final diameters
factor = 1./power
for v, d in diameters.items():
diameters[v] = d**factor
self.error_vertex = error_vertex
if error_vertex:
print('Warnings on diameter for %d vertices'%len(error_vertex))
return diameters
def __interpolate_diameter(self, axe, top_dia, default, order):
top_diameter = self.top_diameter
bottom_diameter = self.bottom_diameter
n = len(axe)
indices = [(i,vid) for i, vid in enumerate(axe) if vid in top_diameter]
# set the default value for the lase node in the axis
last_id = axe[-1]
if last_id not in top_diameter:
klass = self.g.class_name(last_id)
dd = self.dresser.min_topdia.get(klass, default)
indices.append((n-1,last_id))
top_dia[last_id] = dd
if len(indices) > 1:
for i in range(len(indices)-1):
(a, va), (b,vb) = indices[i], indices[i+1]
da, db = top_dia[va], top_dia[vb]
d_step = (db-da)/(b-a)
for index in range(a+1, b):
top_dia[axe[index]] = da + d_step * (index-a)
# Constant interpolation
first_dia = top_dia[indices[0][1]]
for i in range(indices[0][0]):
top_dia[axe[i]] = first_dia
[docs]
def linear_diameter(self, power=2):
""" Traverse all the axes order by order.
The higher order are traversed first.
Retrieve all the defined radius on an axis.
Define a parrameter for each define radius.
Then interpolate and extrapolate
"""
default_diameter = 1 if not self.dresser.min_topdia else min(self.dresser.min_topdia.values())
default_diameter *= 1./self.dresser.diameter_unit
axes = self.axes
top_diameter = self.top_diameter.copy()
if not self.top_diameter and not self.bottom_diameter:
diameter = self.default_algo_diameter(power)
return diameter
orders = list(self.axes)
orders = sorted(orders, reverse=True)
for order in orders:
for axe in self.axes[order]:
self.__interpolate_diameter(axe, top_diameter, default_diameter, order)
return top_diameter
[docs]
def build_mtg_from_radius(self):
""" Decompose the tree (mtg at finest scale) into sub systems
by creating complex which do not have interior nodes with a given values.
Every complex have a defined frontier or are free.
"""
g = self.g
max_scale = g.max_scale()
tree_root = next(g.roots_iter(scale=max_scale))
colors = {}
tv = colors[2] = list(traversal.pre_order(g, tree_root))
bd = self.bottom_diameter
td = self.top_diameter
complex = [v for v in tv if (v == tree_root) or ((v in bd or g.parent(v) in td) and v not in td) ]
colors[1] = complex
mtg, new_map = colored_tree(g, colors)
mtg_root = next(mtg.roots_iter(scale=1))
label = mtg.property('label')
for v in mtg.vertices_iter(scale=1):
label[v] = 'R'+str(v)
return mtg, new_map
[docs]
def decompose_radius(self):
""" Decompose the tree (mtg at finest scale) into sub trees which contains
all the free nodes (nodes without radius).
by creating complex which do not have interior nodes with a given values.
Every tree have a defined frontier or are free.
"""
g = self.g
bd = self.bottom_diameter
td = self.top_diameter
exclude = self.exclude
def exclude_vertex(vid):
return g.class_name(vid) not in exclude
tree_id = {}
trees = {}
max_scale = g.max_scale()
for tree_root in g.roots_iter(scale=max_scale):
for vid in traversal.pre_order2_with_filter(g, tree_root, pre_order_filter=exclude_vertex):
if vid not in td:
pid = g.parent(vid)
if pid in tree_id:
tree = trees[tree_id[pid]]
tree.add_child(pid, vid)
else:
tree = PropertyTree(root=vid)
trees[vid] = tree
tree_id[vid] = tree.root
return trees
#--------------------------------------------------------------------------------------
# algorithms with points
#--------------------------------------------------------------------------------------
#--------------------------------------------------------------------------------------
# algorithms with Length
#--------------------------------------------------------------------------------------
[docs]
def algo_length(self, default_length=None):
"""
Compute the length of each vertex.
The length of a complex is the sum of the length of its components.
Length can be computed from points or is given as a property.
"""
# Check if points are defined
# 1. points not defined
if not self.points:
length = self.algo_length_without_points(default_length=default_length)
else:
length = self._length_from_points()
self._length = length
return self._length
def _length_from_points(self):
g = self.g
origin = self.origin
points = self.points
length = {}
for vid in g.vertices_iter(scale=g.max_scale()):
node = g.node(vid)
# WARNING: Wrong when several trees
pid = g.parent(vid)
pt0 = origin if pid is None else points.get(pid)
pt1 = points.get(vid)
length[vid] = norm(Vector3(*pt1)-Vector3(*pt0))
return length
def _segmentvec_from_points(self):
g = self.g
origin = self.origin
points = self.points
segvec = {}
for vid in g.vertices_iter(scale=g.max_scale()):
node = g.node(vid)
# WARNING: Wrong when several trees
pid = g.parent(vid)
pt0 = origin if pid is None else points.get(pid)
pt1 = points.get(vid)
segvec[vid] = Vector3(*pt1)-Vector3(*pt0)
return segvec
[docs]
def algo_length_without_points(self, default_length=None):
g = self.g
length = self.length
min_length = self.dresser.min_length.copy()
unit = self.dresser.length_unit
if min_length and not default_length:
default_lenght = min(min_length.values())
else:
default_length = default_length if default_length else 1.
_length = {}
_unknows = {}
max_scale = g.max_scale()
# Traverse all the scales.
# in bottom up to propaate values and unknows
for scale in range(max_scale, 1, -1):
for vid in g.vertices_iter(scale=scale):
cid = g.complex(vid)
name = g.class_name(vid)
if vid in length:
_length[cid] = _length.get(cid,0.) + length[vid]
elif vid in _length:
_length[cid] = _length.get(cid,0.) + _length[vid]
elif vid in _unknows:
d = _unknows[vid]
complex_unknow = _unknows.setdefault(cid,{})
for name, x in d.items():
complex_unknow[name] = d.get(name) + x
else:
if name not in min_length:
name = 'default'
_unknows[cid][name] = _unknows.setdefault(cid,{name:0}).get(name,0) + 1
#print _length, _unknows
# Solve the equations in top down.
marked = {}
class Visitor:
def pre_order(self, vid):
if vid in marked or vid in length:
return False
else:
return True
def post_order(self, vid):
pass
visitor = Visitor()
# A possible alternative is to add a weight depending on the ClassName.
for vid in list(length.keys()):
if vid not in _unknows:
continue
n = sum(_unknows.get(vid).values())
x = (length[vid] - _length.get(vid,0)) / float(n)
for v in traversal.pre_order_in_scale(g,vid, visitor):
if v not in length:
print(v)
length[v] = _length.get(v,0) + x * sum(_unknows.get(vid).values())
# Solve all the others without length
for root in g.roots_iter(scale=1):
for v in traversal.pre_order_in_scale(g, root, visitor):
if v not in length:
l = _length.get(v,0)
for name, weight in _unknows.get(v,{}).items():
l+= weight * min_length.get(name, default_length) / unit
length[v] = l
return length
#--------------------------------------------------------------------------------------
# Phyllotaxy algorithms
#--------------------------------------------------------------------------------------
#ToDo
#--------------------------------------------------------------------------------------
# Computing geometric properties
#--------------------------------------------------------------------------------------
[docs]
def compute_length(self,vid= None, default_length = None):
""" Return the length of a vertex.
If no vertex has been defined, return the length of all vertices.
"""
if not self._length:
self.algo_length(default_length)
p = self._length
return p.get(vid) if vid else p
[docs]
def compute_diameter(self, vid=None):
""" Return the diameter of a vertex or for all vertices. """
if not self._diameter:
# get the length and the radius
self.algo_diameter()
p = self._diameter
return p.get(vid) if vid else p
[docs]
def compute_surface(self, vid=None):
""" Return the surface of a vertex or for all vertices.
The surface is for a tapered is
.. math::
\pi (R_{t} + R_b)*h
"""
if not self._surface:
# get the length and the radius
self._surface = surface = {}
diameter = self.compute_diameter()
length = self.compute_length()
g = self.g
for vid in self.g:
if vid not in diameter or vid not in length:
continue
pid = g.parent(vid)
diam_top = diameter[vid]
diam_base = diameter.get(pid,diam_top)
h = length[vid]
s = h*pi*fabs(diam_base+diam_top)/2.
surface[vid] = s
p = self._surface
return p.get(vid) if vid else p
[docs]
def compute_volume(self, vid=None):
""" Return the volume of a vertex or for all vertices.
The volume is for a tapered is
.. math::
V = \frac{\pi h}{3}(R_1^2+R_2^2+R_1 R_2)
"""
if not self._volume:
# get the length and the radius
self._volume = volume= {}
diameter = self.compute_diameter()
length = self.compute_length()
g = self.g
for vid in self.g:
if vid not in diameter or vid not in length:
continue
pid = g.parent(vid)
r_top = diameter[vid]
r_base = diameter.get(pid,r_top)
r_base /= 2.
r_top /= 2.
h = length[vid]
v = (pi*h/3.)*(r_base**2+r_top**2+r_base*r_top)
volume[vid] = v
p = self._volume
return p.get(vid) if vid else p
[docs]
def compute_segmentvec(self, vid=None):
""" Return the segment (i.e. direction with a length) of a vertex or for all vertices. """
if not self._segmentvec:
if self.points:
self._segmentvec = self._segmentvec_from_points()
p = self._segmentvec
return p.get(vid) if vid else p
[docs]
def plot(self, *args, **kwds):
""" Plot a MTG.
:Optional Parameters:
- origins : list of 3D points
- visitor : a function f(g, v, turtle)
This function is called for each vertex of the MTG.
- gc : (bool) generalised cylinder
- turtle: specify the turtle object
"""
g = self.g
# computed properties
diameters = self.compute_diameter()
points = self.points
max_scale = g.max_scale()
origins = kwds.get('origins', [self.origin])
def plantframe_visitor(g, v, turtle):
radius = diameters.get(v)
if radius:
radius = radius /2.
elif g.scale(v) < max_scale:
try:
vm = next(g.component_roots_at_scale_iter(v, scale=max_scale))
radius = diameters.get(vm, 0.)/2.
except:
radius = 0.
if points:
pt = points.get(v)
if not pt and g.scale(v) != max_scale:
try:
vm = next(g.component_roots_at_scale_iter(v, scale=max_scale))
pt = points.get(vm)
except: pass
if pt:
turtle.setId(v)
turtle.lineTo(pt, radius)
else:
turtle.setId(v)
if g.edge_type(v) == '+':
turtle.down()
turtle.F()
if radius:
turtle.setWidth(radius)
turtle.rollL()
visitor = kwds.get('visitor', plantframe_visitor)
gc = kwds.get('gc', True)
_turtle = kwds.get('turtle', None)
scale = kwds.get('scale', max_scale)
if not _turtle:
_turtle = PglTurtle()
roots = g.roots(scale=scale)
for i, rid in enumerate(roots):
if len(origins) > i:
origin = origins[i]
else:
origin = (0,0,0)
_turtle.move(origin)
d = diameters.get(rid)
if not d and scale != max_scale:
try:
vm = next(g.component_roots_at_scale_iter(rid, scale=max_scale))
d = diameters.get(vm)
except:
pass
if d:
_turtle.setWidth(d)
turtle.traverse_with_turtle(g, rid, visitor=visitor, turtle=_turtle, gc=gc)
scene = _turtle.getScene()
if kwds.get('display',True):
Viewer.display(scene)
return scene
[docs]
def plot_property(self, prop, **kwds):
"""
Plot properties of MTG
:Example:
>>> pf.plot_property('length')
"""
pf = self
import matplotlib
import matplotlib.pyplot
props = getattr(self, prop)
pylab_colors = list(matplotlib.colors.cnames.keys())
color = {}
orders = algo.orders(self.g)
color = dict([(k, pylab_colors[orders[k]]) for k in props.keys()])
heights = algo.heights(self.g)
h = dict([(v, heights[v]) for v in props])
for v in props:
matplotlib.pyplot.plot(h[v], props[v], 'o', color = color[v])
#for v in props:
# matplotlib.pyplot.plot(algo.height(pf.g, v), props[v], 'o')
props = getattr(pf, "_"+prop)
h = dict([(v, heights[v]) for v in props])
for order in pf.axes:
for seq in pf.axes[order]:
try:
if order > min(pf.axes.keys()):
x = [h[algo.father(self.g, seq[0])]]
x.extend([h[v] for v in seq])
y = [props[algo.father(self.g, seq[0])]]
y.extend([props[v] for v in seq])
matplotlib.pyplot.plot(x, y)
else:
matplotlib.pyplot.plot([h[v] for v in seq], [props[v] for v in seq])
except:
pass
#matplotlib.pyplot.show()
#return props
###########################################################################################
def iter_order(g, v, edge_type = None):
''' Iter on a tree by considering first
all the vertices of the axe at the first order,
then the vertices at the second order and so on.
'''
pass
def compute_axes(g, v, fixed_points, origin):
marked = {}
axes = {}
others = {}
for vid in traversal.post_order(g,v):
if vid in marked:
continue
_axe = list(simple_axe(g,vid, marked, fixed_points))
_axe.reverse()
_axe, other = list(zip(*_axe))
axes.setdefault(g.order(_axe[0]),[]).append(list(_axe))
others.setdefault(g.order(_axe[0]),[]).append(list(other))
orders = list(axes.keys()).sort()
for order in axes:
for i, axe in enumerate(axes[order]):
other = others[order][i]
_axe = list(zip(axe, other))
if order == 0:
new_points = compute_missing_points(g, _axe, fixed_points, origin=origin)
else:
new_points = compute_missing_points(g, _axe, fixed_points)
fixed_points.update(new_points)
return axes
def simple_axe(g, v, marked, fixed_points):
edge_type = g.property('edge_type')
while v is not None:
if v in fixed_points:
yield v, True
else:
yield v, False
assert v not in marked
marked[v] = True
if g.parent(v) is None or edge_type[v] == '+':
break
v = g.parent(v)
def compute_missing_points(g, axe, fixed_points, origin=None):
"""
axe is a list of tuple containing both vid and a boolean indicated if
a point is defined on a vertex.
return a dict of {id:point}.
"""
new_points = {} # list of points
interval = [] # list of the different interval with missing values
first = False
last = False
prev = -1 # previous id with a defined point value
current = [] # list of points to define an interval
begin = True
for vid, defined in axe:
if begin:
first = not defined
begin = False
if not defined:
if (not current) and (prev != -1) :
current.append(prev)
current.append(vid)
else:
prev = vid
if current:
current.append(vid)
interval.append(current)
current = []
if current:
interval.append(current)
last = True
i0 = []
i1 = []
if first:
i0 = interval[0]
del interval[0]
if last and interval:
i1 = interval[-1]
del interval[-1]
for inter in interval:
n = len(inter)-1
v1, v2 = inter[0], inter[-1]
p1, p2 = fixed_points[v1], fixed_points[v2]
pt12 = Vector3(*p2) -Vector3(*p1)
step = pt12 / n
for i, v in enumerate(inter):
if i == 0:
pt = Vector3(*p1)
elif i == n:
continue
else:
pt = pt + step
new_points[v] = pt
# Management of first and last interval.
# It is a silly algorithm.
if i0 and (i0[-1], True) in axe:
if g.parent(i0[0]) in fixed_points:
origin = fixed_points[g.parent(i0[0])]
if origin :
inter = i0
# The first point do not belong to the interval
n = len(inter)
v2 = inter[-1]
p1, p2 = Vector3(*origin), fixed_points[v2]
pt12 = Vector3(*p2) -Vector3(*p1)
step = pt12 / n
pt = p1
for i, v in enumerate(inter):
if i+1 != n:
pt = pt + step
new_points[v] = pt
else:
"""
# We must compute the min disance between the axe and its parent
v0 = i0[-1]
if (v0, True) in axe:
index0 = axe.index((v0,True))
if index0+1 < len(axe):
v1, defined = axe[index0+1]
if defined:
p2 = Vector3(*fixed_points[v1])
else:
p2 = new_points.get(v1)
if p2:
p1 = Vector3(*fixed_points[v0])
pt21 = Vector3(*p2) -Vector3(*p1)
n = len(i0)-1
for i, v in enumerate(i0[:-1]):
index = n-i
pt = p1 + pt21*(index/n)
#new_points[v] = pt
if v == 2503:
print 2, v, p1, p2, step, pt21
if g.parent(v) in fixed_points:
print step - (p1-fixed_points[g.parent(v)])
"""
print('TODO : point for ', i0[0])
if i1 and (i1 != i0):
v1 = i1[0]
v0 = g.parent(v1)
p1 = fixed_points.get(v0,new_points.get(v0))
if p1 :
p1 = Vector3(*p1)
p2 = fixed_points[v1]
p2 = Vector3(*p2)
pt12 = p2 -p1
n = len(i1)-1
for i, v in enumerate(i1[1:]):
pt = p2 + pt12*(i+1)
new_points[v] = pt
return new_points
def compute_radius(g, v, last_radius):
all_r2 = {}
for vid in traversal.post_order(g, v):
r2 = max(sum([all_r2[c] for c in g.children(vid)]), last_radius)
all_r2[vid] = r2
for k, v in all_r2.items():
all_r2[k] = sqrt(v)
return all_r2
def compute_diameter(g, v, radius, default_value):
all_r = {}
unknow= []
edge_type = g.property('edge_type')
for vid in traversal.post_order(g, v):
if vid in radius:
if radius[vid] < default_value:
all_r[vid] = default_value
else:
all_r[vid] = radius[vid]
elif g.is_leaf(vid):
v = g.complex(vid)
while v:
if v in radius:
all_r[vid] = radius[v]
break
else:
v = g.complex(v)
else:
all_r[vid] = default_value
else:
# pipe model (r_parent **n = sum r_child**n with r ==2)
all_r[vid] = sqrt(sum([all_r[c]**2 for c in g.children(vid)]))
return all_r
def compute_diameter2(g, diameters, property_name='radius',
default_values=None, exclude_symbols=None):
""" Compute the default diameter values of an MTG.
In this algorithm, the MTG is splited into axes. If there are two values for an exis,
these values are interpolated and extrapolated. The interpolation may depends on the
relative position of each internode in the MTG or the length between the internodes.
The default value allows to consider the top of the axis. At the base of the axis,
the extrapolation is linear.
:Properties:
- `g` : the MTG data structure
- `diameters`: a dict containing the diameters for the MTG vertices intially
- `radius`: 3 point to indicate the position of the base of the Tree
- `points`: a dict containing all the 3D coordiantes of the MTG vertices
- `default_radius`: default radius value when there is no diameter given for a given vertex
- `option`: draw all the axes as extrusion shapes or using cylinder
- `colors`: a function returning (r, g, b) values in [0,255]
"""
def build_scene(g, origin, axes, points, diameters, default_radius, option='axe', colors = None, hide=None):
""" Build a scene from the `MTG` g and properties.
Build a 3D scene from the MTG and several properties that have been computed or extracted from the MTG.
Return the 3D PlantGL scene.
:Properties:
- `g` : MTG data structure
- `origin`: 3 point to indicate the position of the base of the Tree
- `points`: a dict containing all the 3D coordiantes of the MTG vertices
- `diameters`: a dict containing the iameters for the MTG vertices
- `default_radius`: default radius value when there is no diameter given for a given vertex
- `option`: draw all the axes as extrusion shapes or using cylinder
- `colors`: a function returning (r, g, b) values in [0,255]
"""
if colors is None:
colors = lambda x: (255, 255, 0)
scene = Scene()
section = Polyline2D.Circle(0.5,10)
polylines = []
radius_law = []
scale = g.max_scale()
rad = diameters
if option == 'cylinder':
for vid in points.keys() :
if hide is not None and hide(vid):
continue
if g.scale(vid) != scale:
continue
parent = g.parent(vid)
if parent is None:
point_parent = origin
elif parent in points:
point_parent = points[parent]
poly = [point_parent, points[vid]]
curve = Polyline(poly)
if vid not in rad or parent not in rad:
rad_vid = rad.get(vid, default_radius)
rad_parent = rad.get(parent, rad_vid) if g.edge_type(vid) == '<' else rad_vid
radius = [[rad_parent]*2, [rad_vid]*2]
color = colors(vid)
if color is None:
color = (255,0,0)
shape = Shape(Extrusion(curve, section, radius), Material(Color3(*color)))
shape.id = vid
scene += shape
else:
rad_vid = rad.get(vid, default_radius)
rad_parent = rad.get(parent, rad_vid) if g.edge_type(vid) == '<' else rad_vid
radius = [[rad_parent]*2, [rad_vid]*2]
color = colors(vid)
if color:
shape = Shape(Extrusion(curve, section, radius),Material(Color3(*color)))
else:
shape = Shape(Extrusion(curve, section, radius))
shape.id = vid
scene += shape
return scene
for order, oaxes in axes.items():
for axe in oaxes:
print(axe)
if hide is not None and hide(axe[0]):
continue
parent = g.parent(axe[0])
if order > 0 and parent and (parent in points):
axe.insert(0,parent)
elif order == 0:
axe.insert(0, axe[0])
poly = [points[vid] for vid in axe if vid in points]
if order == 0:
poly[0] = origin
color = None
for vid in axe:
color = colors(vid)
if color is not None:
break
if not color:
color = (0,0,0)
# Delete null segments
eps = 1
curve = Polyline(poly)
radius = [[rad.get(vid, default_radius)]*2 for vid in axe if vid in points]
if len(radius)>1:
radius[0] = radius[1]
curve, radius = clean_curve(curve, radius)
shape = Shape(Extrusion(curve, section, radius),Material(Color3(*color)))
shape.id = axe[0]
scene += shape
return scene
def debug(g, scene, points, order, scale=3, color=(255,0,0)):
c = Material(Color3(*color))
sphere = Sphere(radius=30)
for id, pt in points.items():
if g.order(id) == order and g.scale(id) == scale:
scene+= Shape(Translated(pt,sphere), c)
return scene
def clean_curve(poly, radius):
""" Remove too small elements.
"""
pts = poly.pointList
n = len(pts)
length = poly.getLength()
mean = length / (n-1)
eps = mean/100.
eps2 = eps**2
index = [True]*n
new_poly = [pts[0]]
new_radius = [radius[0]]
for i in range(1,n):
p12= pts[i] - pts[i-1]
if p12.__normSquared__() > eps2:
new_poly.append(pts[i])
new_radius.append(radius[i])
return Polyline(new_poly), new_radius