# -*- coding: utf-8 -*-
"""
Methods on the mtg for adel programs.
Author: Christophe Pradal
License: CeCILL
Date: 08/07/2008
"""
import csv
from openalea.mtg import MTG, fat_mtg
try:
from openalea.mtg.traversal import *
except:
pass
# from openalea.mtg.io import read_lsystem_string
from openalea.mtg.algo import union
import numpy
try:
from openalea.plantgl.all import (
Scene,
Translation,
Vector3,
Geometry,
AxisRotation,
AxisRotated,
Transform4,
BaseOrientation,
Shape,
Material,
Color3,
PglTurtle,
Mesh,
Translated,
)
except:
Material = tuple
Color3 = tuple
pass
import random
from math import pi
# import read_lsystem string & deps from newmtg/io
try:
from openalea.core.logger import get_logger, logging
logger = get_logger("openalea.mtg")
_ch = logging.StreamHandler()
logger.addHandler(_ch)
except:
logger = None
debug = 0
[docs]
def log(*args):
if debug:
if logger:
logger.debug(" ".join(map(str, args)))
else:
print(" ".join(map(str, args)))
import re
[docs]
def get_expr(s, expr):
res = re.search(expr, s)
_str = ""
if res:
_str = s[res.start() : res.end()]
return _str
[docs]
def get_label(s):
name = r"[a-zA-Z0-9]+"
return get_expr(s, name)
[docs]
def get_name(s):
name = r"[a-zA-Z]+"
return get_expr(s, name)
[docs]
def get_index(s):
name = r"[0-9]+"
return get_expr(s, name)
[docs]
def get_args(s):
args = r"\([0-9,-\.\+]+\)"
return get_expr(s, args)
[docs]
def get_float(s):
args = r"[0-9-\+]+"
num = get_expr(s, args)
return float(num)
[docs]
def read_lsystem_string(string, symbol_at_scale, functional_symbol={}, mtg=None):
"""Read a string generated by a lsystem.
:Parameters:
- `string`: The lsystem string representing the axial tree.
- `symbol_at_scale`: A dict containing the scale for each symbol name.
:Optional parameters:
- `functional_symbol`: A dict containing a function for specific symbols.
The args of the function have to be coherent with those in the string.
The return type of the functions have to be a dictionary of properties: dict(name, value)
:Return:
MTG object
"""
import openalea.plantgl.all as pgl
s = string
def transform(turtle, mesh):
x = turtle.getUp()
z = turtle.getHeading()
bo = pgl.BaseOrientation(x, z ^ x)
matrix = pgl.Transform4(bo.getMatrix())
matrix.translate(turtle.getPosition())
mesh = mesh.transform(matrix)
return mesh
# 1. Create the mtg structure.
if mtg is None:
mtg = MTG()
# 2. add some properties to the MTG
mtg.add_property("index")
mtg.add_property("can_label")
mtg.add_property("geometry")
mtg.add_property("tissue_type")
vid = mtg.root # vid of the support tree, i.e. at the finest scale
current_vertex = mtg.root
branching_stack = []
pending_edge = "" # edge type for the next edge to be created
scale = 0
lsys_symbols = ["[", "]", "/", "+", "^", "f"]
modules = list(symbol_at_scale.keys())
symbols = lsys_symbols + modules
index = dict(list(zip(list(symbol_at_scale.keys()), [0] * len(symbol_at_scale))))
is_ramif = False
# 2. Create a PlantGL Turtle...
turtle = pgl.Turtle()
max_scale = max(symbol_at_scale.values())
for edge_type in symbols:
if edge_type != "f":
s = s.replace(edge_type, "\n%s" % edge_type)
else:
s = s.replace("f(", "\nf(")
l = s.split()
try:
plant_name = [s for s in list(symbol_at_scale.keys()) if "plant" in s.lower()][
0
]
except:
ValueError("""Incorrect plant name (should be plant)""")
for node in l:
# Check if node is a module
tag = node[0]
if tag == "[":
branching_stack.append(vid)
turtle.push()
is_ramif = True
elif tag == "]":
vid = branching_stack.pop()
current_vertex = vid
scale = mtg.scale(vid)
turtle.pop()
is_ramif = False
elif tag == "/":
args = get_args(node[1:])
if args:
angle = get_float(args[1:-1])
turtle.rollR(angle)
else:
turtle.rollR()
elif tag == "+":
args = get_args(node[1:])
if args:
angle = get_float(args[1:-1])
turtle.left(angle)
else:
turtle.left()
elif tag == "^":
args = get_args(node[1:])
if args:
angle = get_float(args[1:-1])
turtle.up(angle)
else:
turtle.up()
elif tag == "f" and node[1] == "(":
args = get_args(node[1:])
if args:
length = get_float(args[1:-1])
if length > 0:
turtle.f(length)
else:
turtle.f()
else:
# add new modules to the mtg (i.e. add nodes)
name = get_name(node)
if name not in modules:
print("Unknow element %s" % name)
continue
module_scale = symbol_at_scale[name]
if is_ramif:
edge_type = "+"
else:
edge_type = "<"
log(node, module_scale, edge_type)
if module_scale == scale:
if mtg.scale(vid) == scale:
vid = mtg.add_child(vid, edge_type=edge_type, label=name)
current_vertex = vid
pending_edge = ""
log(
"",
"Cas 1.1",
scale,
"mtg.scale(vid)",
mtg.scale(vid),
"generated vertex",
vid,
)
assert mtg.scale(vid) == module_scale
else:
# add the edge to the current vertex
current_vertex = mtg.add_child(
current_vertex, edge_type=edge_type, label=name
)
log(
"",
"Cas 1.2",
scale,
"mtg.scale(vid)",
mtg.scale(vid),
"generated vertex",
current_vertex,
)
assert mtg.scale(current_vertex) == module_scale
is_ramif = False
elif module_scale > scale:
log("", "Cas 2", scale, "mtg.scale(vid)", mtg.scale(vid))
old_current_vertex = current_vertex
while module_scale > scale:
if mtg.scale(vid) == scale:
assert vid == current_vertex
vid = mtg.add_component(vid)
current_vertex = vid
log(
"",
"",
"Cas 2.1",
scale,
"generate new component",
current_vertex,
)
scale += 1
if module_scale == scale:
assert mtg.scale(current_vertex) == module_scale
mtg.property("label")[current_vertex] = name
break
else:
scale += 1
current_vertex = mtg.add_component(current_vertex)
else:
log(
node,
"add_child(%d, child=%d)"
% (old_current_vertex, current_vertex),
)
mtg.property("label")[current_vertex] = name
if mtg.scale(vid) == scale:
vid = mtg.add_child(
vid, child=current_vertex, edge_type=edge_type
)
is_ramif = False
else:
assert module_scale < scale
while module_scale < scale:
scale -= 1
current_vertex = mtg.complex(current_vertex)
else:
current_vertex = mtg.add_child(
current_vertex, edge_type=edge_type, label=name
)
assert mtg.scale(current_vertex) == module_scale
# MANAGE the properties, the geometry and the indices!!!
index[name] += 1
if name == plant_name:
for k in list(index.keys()):
if k != name:
index[k] = 0
mtg.property("index")[current_vertex] = index[name]
if name in functional_symbol:
features = eval(node, functional_symbol)
geom = features.get("geometry")
canlabel = features.get("label")
ttype = features.get("tissue_type")
if geom:
# get the transformation from the turtle
geom = transform(turtle, geom)
mtg.property("geometry")[current_vertex] = geom
if name == "StemElement":
# parse args to know how the turtle has to move .
args = get_args(node)[1:-1]
list_args = args.split(",")
length = float(list_args[1]) # 2nd arg
if length > 0:
turtle.f(length)
if canlabel:
canlabel.elt_id = index[name]
plant_id = mtg.complex_at_scale(current_vertex, scale=1)
canlabel.plant_id = mtg.property("index")[plant_id]
mtg.property("can_label")[current_vertex] = canlabel
if ttype:
mtg.property("tissue_type")[current_vertex] = ttype
mtg = fat_mtg(mtg)
return mtg
[docs]
def to_aggregation_table(g):
"""
Convert a mtg `g` to an aggregation table.
Returns a multi-line string with one line per triangle.
"""
# 1. classe at each scale
label = g.property("label")
index = g.property("index")
geometry = g.property("geometry")
can_label = g.property("can_label")
tissue_type = g.property("tissue_type")
symbols = dict(list(zip(iter(label.values()), iter(label.keys()))))
for k, v in symbols.items():
symbols[k] = g.scale(v)
l = list(symbols.items())
l.sort(key=lambda x: x[1])
# scales
header = "Plant Axe Metamer StemElement LeafElement Type"
header = header.split()
# 2. iterate on the geometry at the last scale
assert g.max_scale() == 4
root_axe = next(g.roots_iter(scale=2))
root_metamer = next(g.roots_iter(scale=3))
root_elt = next(g.roots_iter(scale=4))
# compute the number of triangles
nb_lines = sum(mesh.indexListSize() for mesh in geometry.values() if mesh)
lines = numpy.zeros((nb_lines, 6), dtype=int)
# compute relative index for metamer in axe and element in metamer
local_index = {}
# determine the number of element
nb_stem_elements = {}
for root_axe in g.roots_iter(scale=2):
for axe_id in pre_order(g, root_axe):
for i, mid in enumerate(g.components_iter(axe_id)):
local_index[mid] = i + 1
for root_metamer in g.roots_iter(scale=3):
for mid in pre_order(g, root_metamer):
stem_index = 0
leaf_index = 0
for eid in g.components_iter(mid):
if "stem" in label[eid].lower():
stem_index += 1
local_index[eid] = stem_index
else:
leaf_index += 1
local_index[eid] = leaf_index
nb_stem_elements[mid] = stem_index
i = 0
for root_elt in g.roots_iter(scale=4):
for vid in pre_order(g, root_elt):
metamer_id = g.complex(vid)
axe_id = g.complex(metamer_id)
plant_id = g.complex(axe_id)
plant_index = index[plant_id]
axe_index = index[axe_id]
metamer_index = local_index[metamer_id]
element_index = local_index[vid]
stem_index = 0
leaf_index = 0
# element
is_stem = "stem" in label[vid].lower()
if is_stem:
stem_index = element_index
else:
leaf_index = element_index
# determine tissue type
ttype = tissue_type[vid]
##lab = can_label[vid]
## if lab.is_soil():
## ttype = 0
## elif lab.is_leaf():
## ttype = 1#lamina
## elif lab.is_stem():
## if nb_stem_elements[metamer_id] == 1:
## ttype = 2 if lab.optical_id <= 2 else 5#2 = sheath, 5 = ear
## else:
## if element_index == 1:
## ttype = 3 if lab.optical_id <= 2 else 4#internode(3) or peduncle(4)
## else:
## ttype = 2 if lab.optical_id <= 2 else 5#sheath or ear
## else:
## ttype = -1 #unknown
indices = (
plant_index,
axe_index,
metamer_index,
stem_index,
leaf_index,
ttype,
)
# nid = g.node(vid)
# nid.tissue_type = ttype
geom = geometry.get(vid)
if geom:
n = geom.indexListSize()
for j in range(n):
lines[i] = indices
i += 1
lines = lines.transpose()
cols = [c.transpose() for c in lines]
table = dict(list(zip(header, cols)))
return table
[docs]
def to_plantgl(
g,
leaf_material=None,
stem_material=None,
soil_material=None,
colors=None,
ambient_only=False,
):
"""
Returns a plantgl scene from an mtg.
"""
if colors is None:
if leaf_material is None:
leaf_material = Material(Color3(0, 180, 0))
if stem_material is None:
stem_material = Material(Color3(0, 130, 0))
if soil_material is None:
soil_material = Material(Color3(170, 85, 0))
colors = g.property("color")
geometries = g.property("geometry")
labels = g.property("can_label")
scene = Scene()
def geom2shape(vid, mesh, scene):
shape = None
if isinstance(mesh, list):
for m in mesh:
geom2shape(vid, m, scene)
return
if mesh is None:
return
if isinstance(mesh, Shape):
shape = mesh
mesh = mesh.geometry
label = labels.get(vid)
if colors:
if ambient_only:
shape = Shape(
mesh,
Material(
ambient=Color3(*colors.get(vid, [0, 0, 0])),
diffuse=0.0,
specular=Color3(0, 0, 0),
),
)
else:
shape = Shape(mesh, Material(Color3(*colors.get(vid, [0, 0, 0]))))
elif not label:
if not shape:
shape = Shape(mesh)
elif label.is_soil():
shape = Shape(mesh, soil_material)
elif label.is_stem() and label.optical_id <= 1:
shape = Shape(mesh, stem_material)
elif label.is_stem() and label.optical_id > 1:
shape = Shape(mesh, soil_material)
elif label.is_leaf() and label.optical_id <= 1:
shape = Shape(mesh, leaf_material)
elif label.is_leaf() and label.optical_id > 1:
shape = Shape(mesh, soil_material)
shape.id = vid
scene.add(shape)
for vid, mesh in geometries.items():
geom2shape(vid, mesh, scene)
return (scene,)
def _line(ind, pts, label):
s = "p 1 %s 3 %s" % (str(label), " ".join("%.6f" % x for i in ind for x in pts[i]))
return s
[docs]
def to_canestra(g):
"""
Return a string representing a canestra file.
"""
geometry = g.property("geometry")
can_label = g.property("can_label")
begin = "# File generated by OpenAlea.Adel program"
lines = [begin]
max_scale = g.max_scale()
for root_elt in g.roots_iter(scale=max_scale):
for vid in pre_order(g, root_elt):
mesh = geometry.get(vid)
if not mesh:
continue
pts = numpy.array(mesh.pointList, ndmin=2)
indices = numpy.array(mesh.indexList, ndmin=2)
label = can_label[vid]
lines.extend([_line(ind, pts, label) for ind in indices])
lines.append("")
return "\n".join(lines)
[docs]
def planter(g, distribution, random_seed=0, azimuths=None):
"""
Transform a set of plants with given point distributions.
:Parameters:
- g: MTG
- distribution: a list of 2D positions
- random_seed: add a random rotation to each plant if the value is positive.In this case, random_seed is used as a seed.
- azimuth : a list defining plant azimuths (in radians). The list is recycled to achieve positioning of all the plants
"""
# assert g.nb_vertices(scale=1) == len(distribution)
geometry = g.property("geometry")
# store the previous plant translation to not copy each time all the geometry.
if "_plant_translation" not in g.properties():
g.add_property("_plant_translation")
translations = g.property("_plant_translation")
if random_seed > 0:
random.seed(random_seed)
def pt2transfo(pt, previous_pt, rotation):
matrix = Translation(pt).getMatrix()
if rotation != 0:
r = AxisRotation((0, 0, 1), rotation).getMatrix()
matrix = matrix * r * Translation(previous_pt).getMatrix()
return Transform4(matrix), pt
def transform_geom(geom, transfo, translation, rotation):
if isinstance(geom, list):
geom = [transform_geom(g, transfo, translation, rotation) for g in geom]
elif isinstance(geom, Mesh):
geom = geom.transform(transfo) if geom else geom
elif isinstance(geom, Geometry):
geom = Translated(translation, AxisRotated((0, 0, 1), rotation, geom))
elif isinstance(geom, Shape):
geom = Shape(
Translated(translation, AxisRotated((0, 0, 1), rotation, geom.geometry))
)
return geom
plants = g.vertices(scale=1)
plants = plants[: len(distribution)]
max_scale = g.max_scale()
for i, root_elt in enumerate(plants):
previous_translation = translations.get(root_elt, (0, 0, 0))
translations[root_elt] = distribution[i]
# displacement = Vector3(distribution[i])-Vector3(previous_translation)
rotation = 0
if azimuths is not None:
rotation = azimuths[i]
elif random_seed > 0:
rotation = random.random() * pi
transfo, translation = pt2transfo(
Vector3(distribution[i]), -Vector3(previous_translation), rotation
)
l = (
vid
for vid in g.vertices_iter(scale=max_scale)
if g.complex_at_scale(vid, scale=1) == root_elt
)
# for vid in g.components_at_scale(root_elt, 4):
for vid in l:
geom = geometry.get(vid)
geom = transform_geom(geom, transfo, translation, rotation)
if geom:
geometry[vid] = geom
can_label = g.property("can_label").get(vid)
if can_label:
can_label.plant_id = i
[docs]
def duplicate(g, n=1):
if g.nb_vertices(scale=1) == n:
return g
g1 = g.sub_mtg(g.root)
for i in range(n - 1):
g1 = union(g1, g)
return g1
[docs]
class CanMTG(MTG):
def __init__(self, functions, axial_string):
MTG.__init__(self)
symbols = {
"newPlant": 1,
"newAxe": 2,
"newMetamer": 3,
"StemElement": 4,
"LeafElement": 4,
}
g = read_lsystem_string(axial_string, symbols, functions, self)
self.symbols = symbols
[docs]
def apply_property(g, pname, function):
"""Apply a function to each values of a MTG property.
Returns this values as a dict (vid, new value).
"""
prop = g.property(pname)
return dict((k, function(v)) for k, v in prop.items())
CanMTG.planter = planter
CanMTG.to_plantgl = to_plantgl
CanMTG.to_canestra = to_canestra
CanMTG.to_aggregation_table = to_aggregation_table
MTG.planter = planter
MTG.to_plantgl = to_plantgl
MTG.to_canestra = to_canestra
MTG.to_aggregation_table = to_aggregation_table
[docs]
def convert(v, undef="NA"):
res = v
try:
res = int(res)
except ValueError:
try:
res = float(res)
except ValueError:
if res == undef:
res = None
return res
[docs]
def properties(d, exclude=[]):
res = {}
for k, v in d.items():
if k in exclude:
continue
v = convert(v)
if v is not None:
res[k] = v
return res
[docs]
def properties_from_dict(d, index, exclude=None):
if exclude is None:
exclude = []
res = {}
for k in d:
if k in exclude:
continue
res[k] = d[k][index]
return res
[docs]
def topological_table_to_mtg(csv_file, epsilon=1e-6):
f = open(csv_file)
l = f.readline()
sniffer = csv.Sniffer()
dialect = sniffer.sniff(l)
# Go from start of file
f.seek(0)
reader = csv.DictReader(f, dialect=dialect)
# Build the MTG
g = MTG()
topology = ["plant", "axe", "numphy"]
# buffers
prev_plant = 0
prev_axe = -1
prev_metamer = -1
vid_plant = -1
vid_axe = -1
vid_metamer = -1
vid_node = -1
# store the metamer vids of the axis 0
metamers = []
nodes = []
edge_type = "<"
for d in reader:
plant, axe, num_metamer = [int(convert(d.get(x), undef=None)) for x in topology]
# print 'adding vertices for plant: %d, axe:%d, metamer: %d'%(plant, axe, num_metamer)
# Add new plant
if plant != prev_plant:
label = "plant" + str(plant)
vid_plant = g.add_component(g.root, edge_type="/", label=label)
vid_axe = -1
vid_metamer = -1
vid_node = -1
prev_axe = -1
# ??
if num_metamer < prev_metamer:
prev_axe = -1
prev_metamer = -1
# Add an axis
if axe != prev_axe:
label = "axe" + str(axe)
if axe == 0:
vid_axe = g.add_component(vid_plant, edge_type="/", label=label)
vid_node = -1
else:
# args['edge_type'] = '+'
edge_type = "+"
vid_axe = g.add_child(vid_axe, edge_type=edge_type, label=label)
# vid, vid_axe = g.add_child_and_complex(metamers[axe], complex=vid_axe, **args)
# Add metamers
args = properties(d, exclude=topology)
assert num_metamer > 0
label = "metamer" + str(num_metamer)
new_axe = True
if axe == 0 and num_metamer == 1:
metamers = []
edge_type = "/"
vid_metamer = g.add_component(vid_axe, edge_type="/", label=label, **args)
elif num_metamer == 1:
# Add the first element connected to previous metamer on the previous axis
edge_type = "+"
vid_metamer, vid_axe = g.add_child_and_complex(
metamers[axe - 1], complex=vid_axe, edge_type="+", label=label, **args
)
vid_node = nodes[axe - 1]
else:
edge_type = "<"
vid_metamer = g.add_child(vid_metamer, label=label, edge_type="<", **args)
new_axe = False
if axe == 0:
metamers.append(vid_metamer)
# Add internode, sheath, and lamina
nb_internodes = 0
nb_sheath = 0
nb_leaf = 0
# Internode
if args["Ev"] > 0.0:
if args["Esen"] > args["Ev"]:
if vid_node == -1:
vid_node = g.add_component(
vid_metamer,
label="Esen",
edge_type="/",
length=args["Ev"],
po=args["Epos"],
diam=args["Ed"],
)
assert edge_type == "/"
else:
vid_node, vid_metamer = g.add_child_and_complex(
vid_node,
complex=vid_metamer,
label="Esen",
edge_type=edge_type,
length=args["Ev"],
po=args["Epos"],
diam=args["Ed"],
)
else:
if vid_node == -1:
vid_node = g.add_component(
vid_metamer,
label="Egreen",
edge_type="/",
length=args["Ev"] - args["Esen"],
po=args["Epo"],
diam=args["Ed"],
)
assert edge_type == "/"
else:
vid_node, vid_metamer = g.add_child_and_complex(
vid_node,
complex=vid_metamer,
label="Egreen",
edge_type=edge_type,
length=args["Ev"] - args["Esen"],
po=args["Epo"],
diam=args["Ed"],
)
vid_node = g.add_child(
vid_node,
label="Esen",
edge_type="<",
length=args["Esen"],
po=args["Epos"],
diam=args["Ed"],
)
else:
if new_axe and args["Gv"] == 0.0:
if vid_node == -1:
vid_node = g.add_component(
vid_metamer,
label="Egreen",
edge_type="/",
length=0.0,
po=args["Epo"],
diam=args["Ed"],
)
else:
vid_node, vid_metamer = g.add_child_and_complex(
vid_node,
complex=vid_metamer,
label="Egreen",
edge_type=edge_type,
length=0.0,
po=args["Epo"],
diam=args["Ed"],
)
# Sheath
if args["Gv"] > 0.0:
if args["Gsen"] > args["Gv"]:
vid_node = g.add_child(
vid_node,
label="Gsen",
edge_type="<",
length=args["Gv"],
po=args["Gpos"],
diam=args["Gd"],
)
else:
vid_node = g.add_child(
vid_node,
label="Ggreen",
edge_type="<",
length=args["Gv"] - args["Gsen"],
po=args["Gpo"],
diam=args["Gd"],
)
vid_node = g.add_child(
vid_node,
label="Gsen",
edge_type="<",
length=args["Gsen"],
po=args["Gpos"],
diam=args["Gd"],
)
if axe == 0:
nodes.append(vid_node)
# Laminae
if args["Lv"] > 0.0:
if args["Lsen"] > args["Lv"]:
l_node = g.add_child(
vid_node,
label="Lsen",
edge_type="+",
length=args["Lv"],
po=args["Lpos"],
Ll=args["Ll"],
Lw=args["Lw"],
LcType=args["LcType"],
LcIndex=args["LcIndex"],
Linc=args["Linc"],
Laz=args["Laz"],
srb=0,
srt=1,
)
else:
l_node = g.add_child(
vid_node,
label="Lgreen",
edge_type="+",
length=args["Lv"] - args["Lsen"],
po=args["Lpo"],
Ll=args["Ll"],
Lw=args["Lw"],
LcType=args["LcType"],
LcIndex=args["LcIndex"],
Linc=args["Linc"],
Laz=args["Laz"],
srb=0,
srt=1 - (args["Lsen"] / args["Lv"]),
)
l_node = g.add_child(
l_node,
label="Lsen",
edge_type="<",
length=args["Lsen"],
po=args["Lpos"],
Ll=args["Ll"],
Lw=args["Lw"],
LcType=args["LcType"],
LcIndex=args["LcIndex"],
Linc=args["Linc"],
Laz=args["Laz"],
srt=1,
srb=1 - (args["Lsen"] / args["Lv"]),
)
prev_plant = plant
prev_axe = axe
prev_metamer = num_metamer
f.close()
return fat_mtg(g)
[docs]
def internode(g, vid_axe, vid_metamer, prev_node, props, epsilon):
if props["Ev"] < epsilon:
return prev_node
[docs]
def mtg_factory(params, sectors=1):
"""Construct a MTG from a dictionary.
The dictionary contains the parameters and a list of elements. Sector is an integer giving the number of LeafElements per Leaf
The keys of params are:
* plant: index of plant
* axe
* numphy
* Lv
* L_shape
* Lsen
* Lw_shape
* LcType
* LcIndex
* Linc
* Laz
* Lpo
* Lpos
* Gl
* Gv
* Gsen
* Gpo
* Gpos
* Gd
* Ginc
* El
* Ev
* Esen
* Ed
* Einc
* Epo
* Epos
:TODO:
* add length and final_length (DONE)
* fix naming convention for Linc: relative inclination (DONE)
* Add a scale to define the tissue types (ear, sheath, laminae, gain)
* diam and final_diam (resp. width)
* function reset length
* function phenology(g, table) -> dynamic parameters (start_thermaltime, end_thermaltime)
* function growthThermaltime(g, tt, dtt): tt=thermaltime du couvert
* function growthThermaltime(g, tt, dtt, stress factor)
* stress factor: offre/demande
* demand = :math:`D=\int_{tt}^{tt+dtt}{S(x)dx}\rho_s+\int_{tt}^{tt+dtt}{V(x)dx}\rho_v`
* offre = :math:`sum{E_abs}\eps_b`
=> ds = ds_predit* stress_factor
* give the area to the leaf model
* update properties
"""
if not _check_adel_parameters(params):
raise ValueError("Adel parameters are invalid")
g = MTG()
topology = ["plant", "axe", "numphy"]
# buffers
prev_plant = 0
prev_axe = -1
prev_metamer = -1
vid_plant = -1
vid_axe = -1
vid_metamer = -1
vid_node = -1
# store the metamer vids of the axis 0
metamers = []
nodes = []
edge_type = "<"
dp = params
nrow = len(params["plant"])
for i in range(nrow):
# plant, axe, num_metamer = [int(convert(d.get(x),undef=None)) for x in topology]
plant = int(dp["plant"][i])
try:
axe = int(dp["axe"][i])
except:
axe = int(dp["ms_insertion"][i])
num_metamer = int(dp["numphy"][i])
# plant, axe, num_metamer = [int(convert(d.get(x),undef=None)) for x in topology]
# print 'plant: %d, axe:%d, nb_metamers: %d'%(plant, axe, num_metamer)
# Add new plant
if plant != prev_plant:
label = "plant" + str(plant)
vid_plant = g.add_component(g.root, label=label, edge_type="/")
vid_axe = -1
vid_metamer = -1
vid_node = -1
prev_axe = -1
if num_metamer < prev_metamer:
prev_axe = -1
prev_metamer = -1
# Add an axis
if axe != prev_axe:
label = "axe" + str(axe)
if axe == 0:
vid_axe = g.add_component(vid_plant, edge_type="/", label=label)
vid_node = -1
else:
# args['edge_type'] = '+'
edge_type = "+"
assert g.parent(vid_plant) is None
new_axe = g.add_child(vid_axe, edge_type=edge_type, label=label)
print(vid_axe, new_axe, vid_plant)
assert g.parent(vid_plant) is None
# vid, vid_axe = g.add_child_and_complex(metamers[axe], complex=vid_axe, **args)
# Add metamers
args = properties_from_dict(dp, i, exclude=topology)
assert num_metamer > 0
label = "metamer" + str(num_metamer)
if axe == 0 and num_metamer == 1:
metamers = []
edge_type = "/"
vid_metamer = g.add_component(vid_axe, edge_type="/", label=label, **args)
elif num_metamer == 1:
# Add the first element connected to previous metamer on the previous axis
edge_type = "+"
vid_metamer = g.add_component(vid_axe, label=label, **args)
vid_metamer = g.add_child(
metamers[axe - 1], child=vid_metamer, edge_type="+"
)
vid_node = nodes[axe - 1]
else:
edge_type = "<"
vid_metamer = g.add_child(vid_metamer, label=label, edge_type="<", **args)
if axe == 0:
metamers.append(vid_metamer)
# Add internode, sheath, and lamina
# Visible Internode
# if args['Ev'] > 0.:
tissue = "internode"
# Several cases:
Egreen = args["Ev"] - args["Esen"] if args["Ev"] - args["Esen"] > 0.0 else 0.0
Esen = args["Ev"] if Egreen == 0.0 else args["Esen"]
# - Green element and perhaps a senescent one
if vid_node == -1:
# first metamer on the axis
vid_node = g.add_component(
vid_metamer,
label="StemElement",
edge_type="/",
length=Egreen,
final_length=args["El"],
po=args["Epo"],
diam=args["Ed"],
tissue=tissue,
inclination=args["Einc"],
)
assert edge_type == "/"
else:
# first element on the metamer which has a parent metamer on the axis
new_node = g.add_component(vid_metamer)
vid_node = g.add_child(
vid_node,
child=new_node,
label="StemElement",
edge_type=edge_type,
length=Egreen,
final_length=args["El"],
po=args["Epo"],
diam=args["Ed"],
tissue=tissue,
inclination=args["Einc"],
)
# Add a senescent element
# if args['Esen'] > 0.:
vid_node = g.add_child(
vid_node,
label="StemElement",
edge_type="<",
length=Esen,
final_length=args["El"],
po=args["Epos"],
diam=args["Ed"],
tissue=tissue,
sen=True,
inclination=0.0,
)
# Sheath
# Same logic that described previously.
Ggreen = args["Gv"] - args["Gsen"] if args["Gv"] - args["Gsen"] > 0.0 else 0.0
Gsen = args["Gv"] if Ggreen == 0.0 else args["Gsen"]
vid_node = g.add_child(
vid_node,
label="StemElement",
edge_type="<",
length=Ggreen,
final_length=args["Gl"],
po=args["Gpo"],
diam=args["Gd"],
tissue=tissue,
inclination=args["Ginc"],
)
# if args['Gsen'] > 0.:
vid_node = g.add_child(
vid_node,
label="StemElement",
edge_type="<",
length=Gsen,
final_length=args["Gl"],
po=args["Gpos"],
diam=args["Gd"],
sen=True,
tissue=tissue,
inclination=0.0,
)
if axe == 0:
nodes.append(vid_node)
# Lamina
# if args['Lv'] > 0.:
l_node = vid_node
tissue = "lamina"
ds = 1.0 / sectors
srlim = 1.0 - (args["Lsen"] / args["Lv"])
st = ds
edge_type = "+"
Laz = args["Laz"]
Lgreen = args["Lv"] - args["Lsen"] if args["Lv"] - args["Lsen"] > 0.0 else 0.0
Lsen = args["Lv"] if Lgreen == 0.0 else args["Lsen"]
for isect in range(sectors):
# create one green and/or one senescent sector
srb_green = st - ds
srt_green = min(st, max(srb_green, srlim))
srb_sen = srt_green
srt_sen = st
l_node = g.add_child(
l_node,
label="LeafElement",
edge_type=edge_type,
length=args["Lv"],
final_length=args["L_shape"],
po=args["Lpo"],
Ll=args["Ll"],
Lw=args["Lw_shape"],
LcType=args["LcType"],
LcIndex=args["LcIndex"],
Laz=Laz,
Linc=args["Linc"],
srb=srb_green,
srt=srt_green,
tissue=tissue,
)
l_node = g.add_child(
l_node,
label="LeafElement",
edge_type="<",
length=args["Lv"],
final_length=args["L_shape"],
po=args["Lpos"],
Ll=args["Ll"],
Lw=args["Lw_shape"],
LcType=args["LcType"],
LcIndex=args["LcIndex"],
Laz=Laz,
Linc=args["Linc"],
srb=srb_sen,
srt=srt_sen,
sen=True,
tissue=tissue,
)
st += ds
edge_type = "<"
prev_plant = plant
prev_axe = axe
prev_metamer = num_metamer
return fat_mtg(g)
def _check_adel_parameters(params):
return True
[docs]
def compute_element(n, symbols):
leaf = symbols.get("LeafElement")
stem = symbols.get("StemElement")
leaf_rank = int(n.complex().index())
optical_species = int(n.po)
final_length = n.final_length
length = n.length
s_base = n.srb
s_top = n.srt
seed = n.LcIndex
# leaf inclination
linc = n.Linc
if n.label.startswith("L"):
radius_max = n.Lw
element = leaf(
optical_species,
final_length,
length,
radius_max,
s_base,
s_top,
leaf_rank,
seed,
linc,
)
else:
diameter_base = (
n.parent().diam if (n.parent() and n.parent().diam > 0.0) else n.diam
)
diameter_top = n.diam
element = stem(optical_species, length, diameter_base, diameter_top)
can_label = element.get("label")
if can_label:
can_label.elt_id = leaf_rank
plant_node = n.complex_at_scale(scale=1)
can_label.plant_id = plant_node.index()
geom = element.get("geometry")
return geom, can_label
[docs]
def mtg_turtle(g, symbols):
"""Compute the geometry on each node of the MTG using Turtle geometry."""
from openalea.mtg import turtle
plants = g.component_roots_at_scale_iter(g.root, scale=1)
nplants = g.nb_vertices(scale=1)
gt = MTG()
def adel_visitor(g, v, turtle):
# 1. retriev the node
n = g.node(v)
angle = float(n.Laz) if n.Laz else 0.0
turtle.rollL(angle)
if g.edge_type(v) == "+":
if not n.label.startswith("L"):
angle = n.Ginc or n.Einc
angle = float(angle) if angle is not None else 0.0
turtle.up(angle)
# 2. Compute the geometric symbol
mesh, can_label = compute_element(n, symbols)
if mesh:
n.geometry = transform(turtle, mesh)
n.can_label = can_label
# 3. Update the turtle
turtle.setId(v)
if n.label.startswith("S"):
turtle.f(n.length)
# Get the azimuth angle
for plant in plants:
gplant = g.sub_mtg(plant)
scene = turtle.TurtleFrame(gplant, visitor=adel_visitor)
gt = union(gplant, gt)
return gt
[docs]
def mtg_turtle_time(g, symbols, time, update_visitor=None):
"""Compute the geometry on each node of the MTG using Turtle geometry.
Update_visitor is a function called on each node in a pre-order (parent before children).
This function allow to update the parameters and state variables of the vertices.
:Example:
>>> def grow(node, time):
"""
g.properties()["geometry"] = {}
g.properties()["_plant_translation"] = {}
max_scale = g.max_scale()
def compute_element(n, symbols, time):
leaf = symbols.get("LeafElement")
stem = symbols.get("StemElement")
leaf_rank = int(n.complex().index())
optical_species = int(n.po)
metamer = n.complex()
# Length computation
if update_visitor:
length = n.length
final_length = metamer.final_length
else:
final_length = n.final_length
try:
length = (
final_length
* (time - metamer.start_tt)
/ (metamer.end_tt - metamer.start_tt)
if metamer.end_tt and time < metamer.end_tt
else n.length
)
except:
length = n.length
if update_visitor and n.label.startswith("L"):
if metamer.final_length is None:
metamer.final_length = n.final_length
metamer.length = n.length
length = metamer.length
prev_length = (
metamer.final_length
* (n.start_tt - metamer.start_tt)
/ (metamer.end_tt - metamer.start_tt)
)
s_base = (metamer.length - prev_length - n.length) / metamer.length
else:
s_base = n.srb
s_top = n.srt
seed = n.LcIndex
# leaf inclination
if update_visitor and n.label.startswith("L"):
linc = metamer.insertion_angle
else:
linc = n.Linc
element = {}
if n.label.startswith("L"):
radius_max = n.Lw
element = leaf(
optical_species,
final_length,
length,
radius_max,
s_base,
s_top,
leaf_rank,
seed,
linc,
)
else:
diameter_base = (
n.parent().diam if (n.parent() and n.parent().diam > 0.0) else n.diam
)
diameter_top = n.diam
element = stem(optical_species, length, diameter_base, diameter_top)
can_label = element.get("label")
if can_label:
can_label.elt_id = leaf_rank
plant_node = n.complex_at_scale(scale=1)
can_label.plant_id = plant_node.index()
geom = element.get("geometry")
return geom, can_label
def adel_visitor(g, v, turtle, time):
# 1. retriev the node
n = g.node(v)
# Update visitor to compute or modified the node parameters
if update_visitor is not None:
update_visitor(n, time)
if "Leaf" in n.label:
metamer = n.complex()
if (n.start_tt <= time < n.end_tt) or (
(time >= metamer.end_tt) and n.edge_type() == "+"
):
angle = float(metamer.Laz) if metamer.Laz else 0.0
turtle.rollL(angle)
else:
if "Leaf" in n.label:
if n.edge_type() == "+":
angle = float(n.Laz) if n.Laz else 0.0
turtle.rollL(angle)
else:
angle = float(n.Laz) if n.Laz else 0.0
turtle.rollL(angle)
if g.edge_type(v) == "+":
angle = n.Ginc or n.Einc
angle = float(angle) if angle is not None else 0.0
# angle = n.inclination
# angle = float(angle) if angle is not None else 0.
turtle.up(angle)
# 2. Compute the geometric symbol
mesh, can_label = compute_element(n, symbols, time)
if mesh:
n.geometry = transform(turtle, mesh)
n.can_label = can_label
# 3. Update the turtle
turtle.setId(v)
m = n.complex()
if update_visitor:
length = n.length
else:
try:
length = (
n.length * (time - m.start_tt) / (m.end_tt - m.start_tt)
if time < m.end_tt
else n.length
)
except:
length = n.length
if ("Leaf" not in n.label) and (length > 0.0):
turtle.F(length)
# Get the azimuth angle
def traverse_with_turtle_time(g, vid, time, visitor=adel_visitor):
turtle = PglTurtle()
def push_turtle(v):
n = g.node(v)
# if 'Leaf' in n.label:
# return False
try:
start_tt = n.complex().start_tt
if start_tt > time:
return False
except:
pass
if g.edge_type(v) == "+":
turtle.push()
return True
def pop_turtle(v):
n = g.node(v)
try:
start_tt = n.complex().start_tt
if start_tt > time:
return False
except:
pass
if g.edge_type(v) == "+":
turtle.pop()
if g.node(vid).complex().start_tt <= time:
visitor(g, vid, turtle, time)
# turtle.push()
plant_id = g.complex_at_scale(vid, scale=1)
for v in pre_order2_with_filter(g, vid, None, push_turtle, pop_turtle):
if v == vid:
continue
# Done for the leaves
if g.node(v).complex().start_tt > time:
print("Do not consider ", v, time)
continue
visitor(g, v, turtle, time)
scene = turtle.getScene()
return g
for plant_id in g.component_roots_at_scale_iter(g.root, scale=max_scale):
g = traverse_with_turtle_time(g, plant_id, time)
return g
[docs]
def thermal_time(
g, phyllochron=110.0, leaf_duration=1.6, stem_duration=1.6, leaf_falling_rate=10
):
"""
Add dynamic properties on the mtg to simulate developpement
leaf_duration is the phyllochronic time for a leaf to develop from tip appearance to collar appearance
stem_duration is the phyllochronic time for a stem to develop
falling_rate (degrees / phyllochron) is the rate at which leaves fall after colar appearance
"""
plants = g.vertices(scale=1)
metamer_scale = g.max_scale() - 1
for plant in plants:
tt = 0
v = next(g.component_roots_at_scale_iter(plant, scale=metamer_scale))
for metamer in pre_order2(g, v):
end_leaf = tt + phyllochron * leaf_duration
nm = g.node(metamer)
nm.start_tt = tt
nm.end_tt = end_leaf
nm.frate = leaf_falling_rate / phyllochron
sectors = [node for node in nm.components() if "Leaf" in node.label]
stems = [node for node in nm.components() if "Stem" in node.label]
nb_stems = len(stems)
stem_tt = end_leaf
dtt = phyllochron * stem_duration / nb_stems
for stem in stems:
stem.start_tt = stem_tt
stem.end_tt = stem_tt + dtt
stem_tt += dtt
nb_sectors = max(1, len(sectors))
sector_tt = end_leaf
dtt = phyllochron * leaf_duration / nb_sectors
for sector in sectors:
sector.start_tt = sector_tt - dtt
sector.end_tt = sector_tt
sector_tt -= dtt
tt += phyllochron
return g
