using VPL
using Distributions, Plots
# Data types
module TreeTypes
import VPL
# Meristem
struct Meristem <: VPL.Node end
# Bud
struct Bud <: VPL.Node end
# Node
struct Node <: VPL.Node end
# BudNode
struct BudNode <: VPL.Node end
# Internode (needs to be mutable to allow for changes over time)
Base.@kwdef mutable struct Internode <: VPL.Node
length::Float64 = 0.10 # Internodes start at 10 cm
end
# Leaf
Base.@kwdef struct Leaf <: VPL.Node
length::Float64 = 0.20 # Leaves are 20 cm long
width::Float64 = 0.1 # Leaves are 10 cm wide
end
# Graph-level variables
Base.@kwdef struct treeparams
growth::Float64 = 0.1
budbreak::Float64 = 0.25
phyllotaxis::Float64 = 140.0
leaf_angle::Float64 = 30.0
branch_angle::Float64 = 45.0
end
end
import .TreeTypes
# Create geometry + color for the internodes
function VPL.feed!(turtle::Turtle, i::TreeTypes.Internode, vars)
# Rotate turtle around the head to implement elliptical phyllotaxis
rh!(turtle, vars.phyllotaxis)
HollowCylinder!(turtle, length = i.length, height = i.length/15, width = i.length/15,
move = true, color = RGB(0.5,0.4,0.0))
return nothing
end
# Create geometry + color for the leaves
function VPL.feed!(turtle::Turtle, l::TreeTypes.Leaf, vars)
# Rotate turtle around the arm for insertion angle
ra!(turtle, -vars.leaf_angle)
# Generate the leaf
Ellipse!(turtle, length = l.length, width = l.width, move = false,
color = RGB(0.2,0.6,0.2))
# Rotate turtle back to original direction
ra!(turtle, vars.leaf_angle)
return nothing
end
# Insertion angle for the bud nodes
function VPL.feed!(turtle::Turtle, b::TreeTypes.BudNode, vars)
# Rotate turtle around the arm for insertion angle
ra!(turtle, -vars.branch_angle)
end
# Rules
meristem_rule = Rule(TreeTypes.Meristem, rhs = mer -> TreeTypes.Node() +
(TreeTypes.Bud(), TreeTypes.Leaf()) +
TreeTypes.Internode() + TreeTypes.Meristem())
function prob_break(bud)
# We move to parent node in the branch where the bud was created
node = parent(bud)
# We count the number of internodes between node and the first Meristem
# moving down the graph
check, steps = hasDescendent(node, condition = n -> data(n) isa TreeTypes.Meristem)
steps = Int(ceil(steps/2)) # Because it will count both the nodes and the internodes
# Compute probability of bud break and determine whether it happens
if check
prob = min(1.0, steps*vars(bud).budbreak)
return rand() < prob
# If there is no meristem, an error happened since the model does not allow for this
else
error("No meristem found in branch")
end
end
branch_rule = Rule(TreeTypes.Bud,
lhs = prob_break,
rhs = bud -> TreeTypes.BudNode() + TreeTypes.Internode() + TreeTypes.Meristem())
function create_tree(origin, growth, budbreak, orientation)
axiom = T(origin) + RH(orientation) + TreeTypes.Internode() + TreeTypes.Meristem()
tree = Graph(axiom = axiom, rules = (meristem_rule, branch_rule),
vars = TreeTypes.treeparams(growth = growth, budbreak = budbreak))
return tree
end
getInternode = Query(TreeTypes.Internode)
function elongate!(tree, query)
for x in apply(tree, query)
x.length = x.length*(1.0 + vars(tree).growth)
end
end
function growth!(tree, query)
elongate!(tree, query)
rewrite!(tree)
end
function simulate(tree, query, nsteps)
new_tree = deepcopy(tree)
for i in 1:nsteps
growth!(new_tree, query)
end
return new_tree
end
origins = [Vec(i,j,0) for i = 1:2.0:20.0, j = 1:2.0:20.0]
orientations = [rand()*360.0 for i = 1:2.0:20.0, j = 1:2.0:20.0]
growths = rand(LogNormal(-2, 0.3), 10, 10)
budbreaks = rand(Beta(2.0, 10), 10, 10)
forest = vec(create_tree.(origins, growths, budbreaks, orientations));
newforest = [simulate(tree, getInternode, 15) for tree in forest];
render(newforest)Using ray tracer for ground cover
We use the previous simple model of the forest
To use the ray tracer to measure ground cover we just need to use Black materials for the leaves, add a soil tile (also with a black material) and use a directional radiation source that is vertical. The ratio between absorbed power and emmitted power in the soil will give you the ground cover.
Since we are going to be using different types of ray tracers, we wil use the turtle message to use the right material in each case.
function choose_material(message)
if message == :original
nothing
elseif message == :ground_cover
Black()
else
error("The argument message must be :original or :ground_cover")
end
end
# Create geometry + color for the internodes
function VPL.feed!(turtle::Turtle, i::TreeTypes.Internode, vars)
# Rotate turtle around the head to implement elliptical phyllotaxis
rh!(turtle, vars.phyllotaxis)
HollowCylinder!(turtle, length = i.length, height = i.length/15, width = i.length/15,
move = true, color = RGB(0.5,0.4,0.0),
material = choose_material(turtle.message))
return nothing
end
# Create geometry + color for the leaves
function VPL.feed!(turtle::Turtle, l::TreeTypes.Leaf, vars)
# Rotate turtle around the arm for insertion angle
ra!(turtle, -vars.leaf_angle)
# Generate the leaf
Ellipse!(turtle, length = l.length, width = l.width, move = false,
color = RGB(0.2,0.6,0.2), material = choose_material(turtle.message))
# Rotate turtle back to original direction
ra!(turtle, vars.leaf_angle)
return nothing
endWe create a scene with the right message and add a soil tile:
scene = Scene(newforest, message = :ground_cover)
soil = Rectangle(length = 21.0, width = 21.0)
rotatey!(soil, pi/2)
VPL.translate!(soil, Vec(0.0, 10.5, 0.0))
soil_mat = Black()
VPL.add!(scene, mesh = soil, material = soil_mat)We now add a single vertical source on top of the scene
source = DirectionalSource(scene, θ = 0.0, Φ = 0.0, radiosity = 1.0, nrays = 1_000_000)And we run the ray tracer without grid cloner (no need for vertical source and black materials, just ignore the warning!)
settings = RTSettings(nx = 0, ny = 0, parallel = true)
rtobj = RayTracer(scene, source, settings = settings, acceleration = BVH,
rule = SAH{3}(5, 10));
trace!(rtobj)And the ground cover is just the total power absorbed by the soil tile scaled by the emmitted power:
1 - power(soil_mat)[1]/1_000_000/source.power[1]