Skip to content

Commit

Permalink
Merge pull request #3 from VirtualPlantLab/new_mesh
Browse files Browse the repository at this point in the history
Update to new mesh API
  • Loading branch information
AleMorales authored Dec 17, 2024
2 parents 4310971 + 2ca5e00 commit dcb6c3e
Show file tree
Hide file tree
Showing 6 changed files with 71 additions and 70 deletions.
Binary file removed nice_trees.png
Binary file not shown.
32 changes: 16 additions & 16 deletions test/forest.jl
Original file line number Diff line number Diff line change
Expand Up @@ -207,15 +207,15 @@ And we can render the forest with the function `render` as in the binary tree
example but passing the whole forest at once
=#
render(Scene(newforest))
render(Mesh(newforest))
#=
If we iterate 4 more iterations we will start seeing the different individuals
diverging in size due to the differences in growth rates
=#
newforest = [simulate(tree, getInternode, 15) for tree in newforest];
render(Scene(newforest))
render(Mesh(newforest))
#=
## Multithreaded simulation
Expand All @@ -233,7 +233,7 @@ newforest = deepcopy(forest)
@threads for i in eachindex(forest)
newforest[i] = simulate(forest[i], getInternode, 6)
end
render(Scene(newforest, parallel = true))
render(Mesh(newforest, parallel = true))
#=
An alternative way to perform the simulation is to have an outer loop for each timestep and an internal loop over the different trees. Although this approach is not required for this simple model, most FSP models will probably need such a scheme as growth of each individual plant will depend on competition for resources with neighbouring plants. In this case, this approach would look as follows:
Expand All @@ -245,18 +245,18 @@ for step in 1:15
newforest[i] = simulate(newforest[i], getInternode, 1)
end
end
render(Scene(newforest, parallel = true))
render(Mesh(newforest, parallel = true))
#=
# Customizing the scene
# Customizing the mesh
Here we are going to customize the scene of our simulation by adding a horizontal tile represting soil and
Here we are going to customize the mesh of our simulation by adding a horizontal tile represting soil and
tweaking the 3D representation. When we want to combine plants generated from graphs with any other
geometric element it is best to combine all these geometries in a `GLScene` object. We can start the scene
geometric element it is best to combine all these geometries in a `GLScene` object. We can start the mesh
with the `newforest` generated in the above:
=#
scene = Scene(newforest);
mesh = Mesh(newforest);
#=
We can create the soil tile directly, without having to create a graph. The simplest approach is two use
Expand All @@ -273,29 +273,29 @@ rotatey!(soil, pi/2)
VirtualPlantLab.translate!(soil, Vec(0.0, 10.5, 0.0))
#=
We can now add the `soil` to the `scene` object with the `add!` function.
We can now add the `soil` to the `mesh` object with the `add!` function.
=#
VirtualPlantLab.add!(scene, mesh = soil, colors = RGB(1,1,0))
VirtualPlantLab.add!(mesh, soil, colors = RGB(1,1,0))
#=
We can now render the scene that combines the random forest of binary trees and a yellow soil. Notice that
We can now render the mesh that combines the random forest of binary trees and a yellow soil. Notice that
in all previous figures, a coordinate system with grids was being depicted. This is helpful for debugging
your code but also to help setup the scene (e.g. if you are not sure how big the soil tile should be).
your code but also to help setup the mesh (e.g. if you are not sure how big the soil tile should be).
Howver, it may be distracting for the visualization. It turns out that we can turn that off with
`show_axes = false`:
=#
render(scene, axes = false)
render(mesh, axes = false)
#=
We may also want to save a screenshot of the scene. For this, we need to store the output of the `render` function.
We can then resize the window rendering the scene, move around, zoom, etc. When we have a perspective that we like,
We may also want to save a screenshot of the mesh. For this, we need to store the output of the `render` function.
We can then resize the window rendering the mesh, move around, zoom, etc. When we have a perspective that we like,
we can run the `save_scene` function on the object returned from `render`. The argument `resolution` can be adjusted in both
`render` to increase the number of pixels in the final image. A helper function `calculate_resolution` is provided to
compute the resolution from a physical width and height in cm and a dpi (e.g., useful for publications and posters):
=#
res = calculate_resolution(width = 16.0, height = 16.0, dpi = 1_000)
output = render(scene, axes = false, size = res)
output = render(mesh, axes = false, size = res)
export_scene(scene = output, filename = "nice_trees.png")
18 changes: 9 additions & 9 deletions test/growthforest.jl
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ Centre for Crop Systems Analysis - Wageningen University
> - Growth rules, based on information stored in organs (dimensions, carbon assimilation)
> - Update dimensions in function of assimilation
> - Compute sink strength
> - Merge Scenes
> - Merge meshes
> - Generate forest on grid and retrieve canopy-level data (e.g., LAI)
>
Expand Down Expand Up @@ -395,8 +395,8 @@ end
As in the previous example, it makes sense to visualize the forest with a soil
tile beneath it. Unlike in the previous example, we will construct the soil tile
using a dedicated graph and generate a `Scene` object which can later be
merged with the rest of scene generated in daily step:
using a dedicated graph and generate a `Mesh` object which can later be
merged with the rest of mesh generated in daily step:
=#
Base.@kwdef struct Soil <: VirtualPlantLab.Node
Expand All @@ -408,19 +408,19 @@ function VirtualPlantLab.feed!(turtle::Turtle, s::Soil, vars)
end
soil_graph = RA(-90.0) + T(Vec(0.0, 10.0, 0.0)) + ## Moves into position
Soil(length = 20.0, width = 20.0) ## Draws the soil tile
soil = Scene(Graph(axiom = soil_graph));
soil = Mesh(Graph(axiom = soil_graph));
render(soil, axes = false)
#=
And the following function renders the entire scene (notice that we need to
use `display()` to force the rendering of the scene when called within a loop
And the following function renders the entire mesh (notice that we need to
use `display()` to force the rendering of the mesh when called within a loop
or a function):
=#
function render_forest(forest, soil)
scene = Scene(vec(forest)) ## create scene from forest
scene = Scene([scene, soil]) ## merges the two scenes
render(scene)
mesh = Mesh(vec(forest)) ## create mesh from forest
mesh = Mesh([mesh, soil]) ## merges the two scenes
render(mesh)
end
#=
Expand Down
65 changes: 33 additions & 32 deletions test/raytracedforest.jl
Original file line number Diff line number Diff line change
Expand Up @@ -206,14 +206,14 @@ end
As growth is now dependent on intercepted PAR via RUE, we now need to simulate
light interception by the trees. We will use a ray-tracing approach to do so.
The first step is to create a scene with the trees and the light sources. As for
rendering, the scene can be created from the `forest` object by simply calling
`Scene(forest)` that will generate the 3D meshes and connect them to their
The first step is to create a mesh with the trees and the light sources. As for
rendering, the mesh can be created from the `forest` object by simply calling
`Mesh(forest)` that will generate the 3D meshes and connect them to their
optical properties.
However, we also want to add the soil surface as this will affect the light
distribution within the scene due to reflection from the soil surface. This is
similar to the customized scene that we created before for rendering, but now
distribution within the mesh due to reflection from the soil surface. This is
similar to the customized mesh that we created before for rendering, but now
for the light simulation.
=#
Expand All @@ -225,29 +225,29 @@ function create_soil()
end
function create_scene(forest)
# These are the trees
scene = Scene(vec(forest))
mesh = Mesh(vec(forest))
# Add a soil surface
soil = create_soil()
soil_material = Lambertian= 0.0, ρ = 0.21)
add!(scene, mesh = soil, materials = soil_material)
# Return the scene
return scene
add!(mesh, soil, materials = soil_material)
# Return the mesh
return mesh
end
#=
Given the scene, we can create the light sources that can approximate the solar
Given the mesh, we can create the light sources that can approximate the solar
irradiance on a given day, location and time of the day using the functions from
the package (see package documentation for details). Given the latitude,
day of year and fraction of the day (`f = 0` being sunrise and `f = 1` being sunset),
the function `clear_sky()` computes the direct and diffuse solar radiation assuming
a clear sky. These values may be converted to different wavebands and units using
`waveband_conversion()`. Finally, the collection of light sources approximating
the solar irradiance distribution over the sky hemisphere is constructed with the
function `sky()` (this last step requires the 3D scene as input in order to place
function `sky()` (this last step requires the 3D mesh as input in order to place
the light sources adequately).
=#
function create_sky(;scene, lat = 52.0*π/180.0, DOY = 182)
function create_sky(;mesh, lat = 52.0*π/180.0, DOY = 182)
# Fraction of the day and day length
fs = collect(0.1:0.1:0.9)
dec = declination(DOY)
Expand All @@ -264,7 +264,7 @@ function create_sky(;scene, lat = 52.0*π/180.0, DOY = 182)
Idir_PAR = f_dir.*Idir
Idif_PAR = f_dif.*Idif
# Create the dome of diffuse light
dome = sky(scene,
dome = sky(mesh,
Idir = 0.0, ## No direct solar radiation
Idif = sum(Idir_PAR)/10*DL, ## Daily Diffuse solar radiation
nrays_dif = 1_000_000, ## Total number of rays for diffuse solar radiation
Expand All @@ -274,20 +274,20 @@ function create_sky(;scene, lat = 52.0*π/180.0, DOY = 182)
nphi = 12) ## Number of discretization steps in the azimuth angle
# Add direct sources for different times of the day
for I in Idir_PAR
push!(dome, sky(scene, Idir = I/10*DL, nrays_dir = 100_000, Idif = 0.0)[1])
push!(dome, sky(mesh, Idir = I/10*DL, nrays_dir = 100_000, Idif = 0.0)[1])
end
return dome
end
#=
The 3D scene and the light sources are then combined into a `RayTracer` object,
The 3D mesh and the light sources are then combined into a `RayTracer` object,
together with general settings for the ray tracing simulation chosen via `RTSettings()`.
The most important settings refer to the Russian roulette system and the grid
cloner (see section on Ray Tracing for details). The settings for the Russian
roulette system include the number of times a ray will be traced
deterministically (`maxiter`) and the probability that a ray that exceeds `maxiter`
is terminated (`pkill`). The grid cloner is used to approximate an infinite canopy
by replicating the scene in the different directions (`nx` and `ny` being the
by replicating the mesh in the different directions (`nx` and `ny` being the
number of replicates in each direction along the x and y axes, respectively). It
is also possible to turn on parallelization of the ray tracing simulation by
setting `parallel = true` (currently this uses Julia's builtin multithreading
Expand All @@ -296,26 +296,26 @@ capabilities).
In addition `RTSettings()`, an acceleration structure and a splitting rule can
be defined when creating the `RayTracer` object (see ray tracing documentation
for details). The acceleration structure allows speeding up the ray tracing
by avoiding testing all rays against all objects in the scene.
by avoiding testing all rays against all objects in the mesh.
=#
function create_raytracer(scene, sources)
function create_raytracer(mesh, sources)
settings = RTSettings(pkill = 0.9, maxiter = 4, nx = 5, ny = 5, parallel = true)
RayTracer(scene, sources, settings = settings, acceleration = BVH,
RayTracer(mesh, sources, settings = settings, acceleration = BVH,
rule = SAH{3}(5, 10));
end
#=
The actual ray tracing simulation is performed by calling the `trace!()` method
on the ray tracing object. This will trace all rays from all light sources and
update the radiant power absorbed by the different surfaces in the scene inside
update the radiant power absorbed by the different surfaces in the mesh inside
the `Material` objects (see `feed!()` above):
=#
function run_raytracer!(forest; DOY = 182)
scene = create_scene(forest)
sources = create_sky(scene = scene, DOY = DOY)
rtobj = create_raytracer(scene, sources)
mesh = create_scene(forest)
sources = create_sky(mesh = mesh, DOY = DOY)
rtobj = create_raytracer(mesh, sources)
trace!(rtobj)
return nothing
end
Expand Down Expand Up @@ -544,32 +544,33 @@ end
As in the previous example, it makes sense to visualize the forest with a soil
tile beneath it. Unlike in the previous example, we will construct the soil tile
using a dedicated graph and generate a `Scene` object which can later be
merged with the rest of scene generated in daily step:
using a dedicated graph and generate a `Mesh` object which can later be
merged with the rest of mesh generated in daily step:
=#
Base.@kwdef struct Soil <: VirtualPlantLab.Node
length::Float64
width::Float64
end
function VirtualPlantLab.feed!(turtle::Turtle, s::Soil, data)
Rectangle!(turtle, length = s.length, width = s.width, colors = RGB(255/255, 236/255, 179/255))
Rectangle!(turtle, length = s.length, width = s.width, colors = RGB(255/255, 236/255, 179/255),
materials = Lambertian= 0.0, ρ = 0.21))
end
soil_graph = RA(-90.0) + T(Vec(0.0, 10.0, 0.0)) + ## Moves into position
Soil(length = 20.0, width = 20.0) ## Draws the soil tile
soil = Scene(Graph(axiom = soil_graph));
soil = Mesh(Graph(axiom = soil_graph));
render(soil, axes = false)
#=
And the following function renders the entire scene (notice that we need to
use `display()` to force the rendering of the scene when called within a loop
And the following function renders the entire mesh (notice that we need to
use `display()` to force the rendering of the mesh when called within a loop
or a function):
=#
function render_forest(forest, soil)
scene = Scene(vec(forest)) ## create scene from forest
scene = Scene([scene, soil]) ## merges the two scenes
render(scene)
mesh = Mesh(vec(forest)) ## create mesh from forest
mesh = Mesh([mesh, soil]) ## merges the two scenes
render(mesh)
end
#=
Expand Down
22 changes: 11 additions & 11 deletions test/snowflakes.jl
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@ class to implement the edges of the snowflake. This can be achieved as follows:
=#
using VirtualPlantLab
import GLMakie ## Import rather than "using" to avoid masking Scene
import GLMakie ## Import rather than "using" to avoid masking Mesh
using ColorTypes ## To define colors for the rendering
module sn
import VirtualPlantLab
Expand Down Expand Up @@ -128,7 +128,7 @@ After defining the method, we can now call the function render on the graph to
generate a 3D interactive image of the Koch snowflake in the current state
=#
sc = Scene(Koch)
sc = Mesh(Koch)
render(sc, axes = false)
#=
Expand All @@ -138,7 +138,7 @@ snowflake. Let's execute the rules once to verify that we get the 2nd iteration
=#
rewrite!(Koch)
render(Scene(Koch), axes = false)
render(Mesh(Koch), axes = false)
#=
And two more times
Expand All @@ -147,7 +147,7 @@ And two more times
for i in 1:3
rewrite!(Koch)
end
render(Scene(Koch), axes = false)
render(Mesh(Koch), axes = false)
#=
# Other snowflake fractals
Expand All @@ -172,13 +172,13 @@ iterations and render the results
=#
## First iteration
rewrite!(Koch2)
render(Scene(Koch2), axes = false)
render(Mesh(Koch2), axes = false)
## Second iteration
rewrite!(Koch2)
render(Scene(Koch2), axes = false)
render(Mesh(Koch2), axes = false)
## Third iteration
rewrite!(Koch2)
render(Scene(Koch2), axes = false)
render(Mesh(Koch2), axes = false)
#=
This is know as [Koch
Expand All @@ -190,18 +190,18 @@ axiom:
=#
axiomCesaro = sn.E(L) + RU(90.0) + sn.E(L) + RU(90.0) + sn.E(L) + RU(90.0) + sn.E(L)
Cesaro = Graph(axiom = axiomCesaro, rules = (rule2,))
render(Scene(Cesaro), axes = false)
render(Mesh(Cesaro), axes = false)
#=
And, as before, let's go through the first three iterations
=#
## First iteration
rewrite!(Cesaro)
render(Scene(Cesaro), axes = false)
render(Mesh(Cesaro), axes = false)
## Second iteration
rewrite!(Cesaro)
render(Scene(Cesaro), axes = false)
render(Mesh(Cesaro), axes = false)
## Third iteration
rewrite!(Cesaro)
render(Scene(Cesaro), axes = false)
render(Mesh(Cesaro), axes = false)
4 changes: 2 additions & 2 deletions test/tree.jl
Original file line number Diff line number Diff line change
Expand Up @@ -186,11 +186,11 @@ newtree = simulate(tree, getInternode, 2)
The binary tree after two iterations has two branches, as expected:
=#
render(Scene(newtree))
render(Mesh(newtree))
#=
Notice how the lengths of the prisms representing internodes decreases as the branching order increases, as the internodes are younger (i.e. were generated fewer generations ago). Further steps will generate a structure that is more tree-like.
=#
newtree = simulate(newtree, getInternode, 15)
render(Scene(newtree))
render(Mesh(newtree))

0 comments on commit dcb6c3e

Please sign in to comment.