Adapt to Meshes v0.44

Project.toml

@@ -17,6 +17,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 TransformsBase = "28dd2a49-a57a-4bfb-84ca-1a49db9b96b8"
 UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [weakdeps]
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
@@ -28,8 +29,8 @@ PlantGeomMakie = "Makie"
 ColorSchemes = "3"
 Colors = "0.12"
 EzXML = "1"
-Makie = "0.15, 0.16, 0.17, 0.18, 0.19, 0.20"
-Meshes = "0.42"
+Makie = "0.21"
+Meshes = "0.44"
 MultiScaleTreeGraph = "0.14"
 Observables = "0.5"
 RecipesBase = "1"

docs/Project.toml

@@ -2,6 +2,7 @@
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 JSServe = "824d6782-a2ef-11e9-3a09-e5662e0c26f9"
+Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa"
 MultiScaleTreeGraph = "dd4a991b-8a45-4075-bede-262ee62d5583"
 PlantGeom = "5edaa67e-25db-4eb9-bf81-05d793b2238d"
 PlotlyJS = "f0f68f2c-4968-5e81-91da-67840de0976a"

docs/make.jl

@@ -1,4 +1,5 @@
 using CairoMakie
+using Meshes
 using PlantGeom
 using Documenter
 

docs/src/index.md

@@ -17,16 +17,15 @@ The package provides different functionalities, the main ones being:
 Note that PlantGeom reserves the `:geometry` attribute in the nodes (*e.g.* organs). It uses it to store the 3D geometry as a special structure ([`geometry`](@ref)).
 
 ```@setup animation
-using CairoMakie, PlantGeom, MultiScaleTreeGraph # Note: CairoMakie must be loaded before PlantGeom to access the extensions
+using CairoMakie, Meshes, PlantGeom, MultiScaleTreeGraph # Note: CairoMakie must be loaded before PlantGeom to access the extensions
 opf = read_opf(joinpath(dirname(dirname(pathof(PlantGeom))),"test","files","simple_plant.opf"))
 
 # First, we compute the 3D coordinates for each node in the MTG:
 transform!(opf, refmesh_to_mesh!)
 # And compute the max z of each node based on their mesh:
 transform!(opf, zmax => :z_node, ignore_nothing = true)
 # Or the z coordinate of each vertez of each node mesh:
-transform!(opf, :geometry => (x -> [i.coords[3] for i in x.mesh.vertices]) => :z_vertex, ignore_nothing = true)
-
+transform!(opf, :geometry => (x -> [Meshes.coords(i).z for i in Meshes.vertices(x.mesh)]) => :z_vertex, ignore_nothing = true)
 
 # Then we make a Makie figure:
 f = Figure()
@@ -57,15 +56,15 @@ end
 If you want to reproduce the animation, you can look at the code below. Otherwise, please head to the next section.
 
 ```julia
-using PlantGeom, MultiScaleTreeGraph, CairoMakie
+using CairoMakie, Meshes, PlantGeom, MultiScaleTreeGraph
 opf = read_opf(joinpath(dirname(dirname(pathof(PlantGeom))),"test","files","simple_plant.opf"))
 
 # First, we compute the 3D coordinates for each node in the MTG:
 transform!(opf, refmesh_to_mesh!)
 # And compute the max z of each node based on their mesh:
 transform!(opf, zmax => :z_node, ignore_nothing = true)
 # Or the z coordinate of each vertez of each node mesh:
-transform!(opf, :geometry => (x -> [i.coords[3] for i in x.mesh.vertices]) => :z_vertex, ignore_nothing = true)
+transform!(opf, :geometry => (x -> [Meshes.coords(i)z for i in Meshes.vertices(x.mesh)]) => :z_vertex, ignore_nothing = true)
 
 
 # Then we make a Makie figure:

docs/src/makie_3d.md

@@ -167,7 +167,7 @@ f
 
 ```julia
 # Compute the z position of each vertices in each mesh:
-transform!(mtg, :geometry => (x -> [i.coords[3] for i in x.mesh.vertices]) => :z, ignore_nothing = true)
+transform!(mtg, :geometry => (x -> [Meshes.coords(i).z for i in Meshes.vertices(x.mesh)]) => :z, ignore_nothing = true)
 viz(mtg, color = :z, showfacets = true, color_vertex = true)
 ```
 

ext/PlantGeomMakie.jl

@@ -11,6 +11,7 @@ import MultiScaleTreeGraph
 import MultiScaleTreeGraph: get_attributes
 import ColorSchemes: get, rainbow, colorschemes, ColorScheme
 import UUIDs
+import Unitful
 
 MeshesMakieExt = Base.get_extension(Meshes, :MeshesMakieExt)
 

ext/makie_recipes/colorbar.jl

@@ -37,18 +37,19 @@ function PlantGeom.colorbar(parent, plotobject; kwargs...)
     end
 
     # Because we extend the `Viz` type, we need to check if the user has given a color range.
-    # If we defined our own e.g. `PlantViz` type, we could have defined a `color_range` field in it directly.
-    if hasproperty(plotobject.attributes, :color_range)
-        if isa(plotobject.attributes.color_range, Observables.Observable)
-            colorbar_limits = plotobject.attributes.color_range
+    # If we defined our own e.g. `PlantViz` type, we could have defined a `colorrange` field in it directly.
+
+    if hasproperty(plotobject.attributes, :colorrange) && (!isa(plotobject.attributes.colorrange, Observables.Observable) || plotobject.attributes.colorrange[] !== nothing)
+        if isa(plotobject.attributes.colorrange, Observables.Observable)
+            # colorbar_limits = Unitful.ustrip.(plotobject.attributes.colorrange[])
+            colorbar_limits = Makie.lift(x -> Unitful.ustrip.(x), plotobject.attributes.colorrange)
         else
-            colorbar_limits = Observables.Observable(plotobject.attributes.color_range)
+            colorbar_limits = Observables.Observable(plotobject.attributes.colorrange)
         end
     else
         # Get the attribute values without nothing values:    
-        colorbar_limits = Makie.@lift PlantGeom.attribute_range($mtg, $color)
+        colorbar_limits = Makie.@lift PlantGeom.attribute_range($mtg, $color, ustrip=true)
     end
-
     colormap = Makie.lift(get_colormap, plotobject.attributes.colormap)
 
     Makie.Colorbar(

ext/makie_recipes/opf/plot_opf.jl

@@ -12,7 +12,7 @@ The plot object can have the following optional arguments:
 - `color`: The color to be used for the plot. Can be a colorant, an attribute of the MTG, or a dictionary of colors for each reference mesh.
 - `alpha`: The alpha value to be used for the plot. Should be a float between 0 and 1.
 - `colormap`: The colorscheme to be used for the plot. Can be a Symbol or a ColorScheme.
-- `color_range`: The range of values to be used for the colormap. Should be a tuple of floats.
+- `colorrange`: The range of values to be used for the colormap. Should be a tuple of floats (optionally with units if *e.g.* z position).
 - `segmentcolor`: The color to be used for the facets. Should be a colorant or a symbol of color.
 - `showsegments`: A boolean indicating whether the facets should be shown or not.
 - `segmentsize`: The size of the segments. Should be a float.
@@ -129,7 +129,7 @@ function plot_opf(colorant::Observables.Observable{AttributeColorant}, plot)
     # from the plot. Instead, we need to check here if the argument is given, and give the default
     # value if not.
     # Note: If we defined our own e.g. `PlantViz` type, we could have defined a `color_missing` and 
-    # `color_range` fields in it directly.
+    # `colorrange` fields in it directly.
 
     # Are the colors given for each vertex in the meshes, or for each reference mesh?
     # Note that we can have several values if we have several timesteps too.
@@ -145,11 +145,11 @@ function plot_opf(colorant::Observables.Observable{AttributeColorant}, plot)
         # Get the attribute values without nothing values:    
         color_missing = Observables.Observable(RGBA(0, 0, 0, 0.3))
     end
-    if hasproperty(plot, :color_range)
-        color_range = plot[:color_range]
+    if hasproperty(plot, :colorrange) && (!isa(plot[:colorrange], Observables.Observable) || plot[:colorrange][] !== nothing)
+        color_range = plot[:colorrange]
     else
         # Get the attribute values without nothing values:    
-        color_range = Makie.@lift PlantGeom.attribute_range($opf, $colorant)
+        color_range = Makie.@lift PlantGeom.attribute_range($opf, $colorant, ustrip=true)
     end
 
     if hasproperty(plot, :index)

ext/makie_recipes/opf_recipe.jl

@@ -66,7 +66,7 @@ vertex_color = get_color(1:nvertices(get_ref_meshes(opf))[1], [1,nvertices(get_r
 viz(opf, color = Dict(1 => vertex_color))
 
 # Or even coloring by the value of the Z coordinates of each vertex:
-transform!(opf, :geometry => (x -> [i.coords[3] for i in x.mesh.vertices]) => :z, ignore_nothing = true)
+transform!(opf, :geometry => (x -> [Meshes.coords(i).z for i in Meshes.vertices(x.mesh)]) => :z, ignore_nothing = true)
 viz(opf, color = :z, showsegments = true)
 
 f,a,p = viz(opf, color = :z, showsegments = true)

src/PlantGeom.jl

@@ -4,14 +4,15 @@ using MultiScaleTreeGraph
 import Observables # For to_value (get an observable value)
 # For 3D (OPF):
 import Meshes
-import Meshes: Translate, Affine, Rotate, Scale, Vec3
+import Meshes: Translate, Affine, Rotate, Scale, Vec
 import Meshes: viz, viz!
 import TransformsBase: parameters, Identity, Transform, →, SequentialTransform
 import TransformsBase: isrevertible, isinvertible
 import TransformsBase: apply, revert, reapply, inverse
 import TransformsBase: parameters
 import Rotations: Rotation, RotZ
-
+import Unitful
+import Unitful: @u_str
 import Tables
 
 # import GeometryBasics

src/colors/colors.jl

@@ -41,18 +41,20 @@ get_color(1:2, 1:10, colormap = colorschemes[:viridis]) # returns RGB{N0f8}(0.26
 get_color(1:2, 1:10, 1, colormap = colorschemes[:viridis]) # returns RGB{N0f8}(0.267004,0.00487433,0.329415)
 """
 function get_color(var::T, range_var, index::Nothing=nothing; colormap=colorschemes[:viridis]) where {T<:AbstractArray}
+    range_var = Unitful.ustrip.(range_var)
     x2 = (range_var[2] - range_var[1])
     # get the color based on a colormap and the normalized attribute value
-    [get(colormap, (i - range_var[1]) / x2) for i in var]
+    [get(colormap, (i - range_var[1]) / x2) for i in Unitful.ustrip.(var)]
 end
 
 function get_color(var::T, range_var, index::I; colormap=colorschemes[:viridis]) where {T<:AbstractArray,I<:Integer}
     get_color(var[index], range_var; colormap=colormap)
 end
 
 function get_color(var, range_var, index::I=1; colormap=colorschemes[:viridis]) where {I<:Integer}
+    range_var = Unitful.ustrip.(range_var)
     # get the color based on a colormap and the normalized attribute value
-    get(colormap, (var - range_var[1]) / (range_var[2] - range_var[1]))
+    get(colormap, (Unitful.ustrip(var) - range_var[1]) / (range_var[2] - range_var[1]))
 end
 
 

src/equality.jl

@@ -31,7 +31,7 @@ Test RefMeshes equality.
 function Base.:(==)(a::T, b::T) where {T<:geometry}
     isequal(a.ref_mesh, b.ref_mesh) &&
         isequal(a.ref_mesh_index, b.ref_mesh_index) &&
-        isequal(a.transformation(Meshes.Point3(1, 2, 3)), b.transformation(Meshes.Point3(1, 2, 3))) &&
+        isequal(a.transformation(Meshes.Point(1, 2, 3)), b.transformation(Meshes.Point(1, 2, 3))) &&
         # NB: transform a point here because transformations can't be compared directly
         isequal(a.dUp, b.dUp) &&
         isequal(a.dDwn, b.dDwn) &&

src/helpers.jl

@@ -19,7 +19,7 @@ end
 function normals(mesh::RefMesh{S,ME,M,N,T}) where {S,ME<:Meshes.SimpleMesh,M,N,T}
     if length(mesh.normals) == 0
         return SVector{length(Meshes.topology(mesh.mesh).connec)}(
-            Meshes.Point3(Meshes.normal(Meshes.Triangle(mesh.mesh.vertices[[tri.indices...]]))) for tri in Meshes.topology(mesh.mesh).connec
+            Meshes.Point(Meshes.normal(Meshes.Triangle(Meshes.vertices(mesh.mesh)[[tri.indices...]]))) for tri in Meshes.topology(mesh.mesh).connec
         )
     else
         return mesh.normals
@@ -38,7 +38,7 @@ Compute per vertex normals and return them as a `StaticArrays.SVector`.
 # https://stackoverflow.com/questions/45477806/general-method-for-calculating-smooth-vertex-normals-with-100-smoothness?noredirect=1&lq=1
 """
 function normals_vertex(mesh::RefMesh)
-    vertex_normals = fill(Meshes.Point3(0.0, 0.0, 0.0), Meshes.nvertices(mesh))
+    vertex_normals = fill(Meshes.Vec(0.0, 0.0, 0.0), Meshes.nvertices(mesh))
     for (i, tri) in enumerate(Meshes.topology(mesh.mesh).connec)
         vertex_normals[tri.indices[1]] = mesh.normals[i]
         vertex_normals[tri.indices[2]] = mesh.normals[i]
@@ -49,7 +49,7 @@ function normals_vertex(mesh::RefMesh)
 end
 
 function normals_vertex(mesh::Meshes.SimpleMesh)
-    vertex_normals = fill(Meshes.Point3(0.0, 0.0, 0.0), Meshes.nvertices(mesh))
+    vertex_normals = fill(Meshes.Vec(0.0, 0.0, 0.0), Meshes.nvertices(mesh))
     for (i, tri) in enumerate(Meshes.topology(mesh).connec)
         tri_norm = Meshes.normal(Meshes.Triangle(mesh.vertices[[tri.indices...]]...))
         vertex_normals[tri.indices[1]] = tri_norm

src/meshes/summary_coordinates.jl

@@ -37,5 +37,5 @@ Apply function `f` over the mesh coordinates `coord`.
 Values for `coord` can be 1 for x, 2 for y and 3 for z.
 """
 function map_coord(f, mesh, coord)
-    f([i.coords[coord] for i in mesh.vertices])
+    f([Meshes.to(i)[coord] for i in Meshes.vertices(mesh)])
 end

src/meshes/transformations.jl

@@ -51,7 +51,7 @@ leaf_node = get_node(opf4, 8)
 parent_zmax = zmax(leaf_node.parent)
 
 # Define a rotation of the mesh around the Z axis defined by the parent node max Z:
-transformation = recenter(Rotate(RotZ(1.0)), Point3(0.0, 0.0, parent_zmax))
+transformation = recenter(Rotate(RotZ(1.0)), Point(0.0, 0.0, parent_zmax))
 
 # Update the transformation matrix of the leaf and its mesh:
 transform_mesh!(leaf_node, transformation)

src/opf/mtg_recipe_helpers.jl

@@ -164,7 +164,7 @@ function mtg_coordinates_df!(mtg, attr=:YY; force=false)
     DataFrame(mtg, unique([:XX, :YY, :ZZ, :XX_from, :YY_from, :ZZ_from, attr]))
 end
 
-function attribute_range(mtg, attr)
+function attribute_range(mtg, attr; ustrip=false)
     attr_name = attr_colorant_name(attr)
     vals =
         MultiScaleTreeGraph.descendants(
@@ -183,6 +183,10 @@ function attribute_range(mtg, attr)
         range_val = (minimum(minimum.(vals_no_nothing)), maximum(maximum.(vals_no_nothing)))
     end
 
+    if ustrip
+        return Unitful.ustrip.(range_val)
+    end
+
     return range_val
 end
 

src/opf/read_opf.jl

@@ -178,23 +178,31 @@ function parse_meshBDD!(node)
                 content = parse_opf_array(i.content, Int) .+ 1
                 # NB: adding 1 to the faces because the opf is 0-based but Julia is 1-based
 
-                faces3d = Meshes.Connectivity[
+                faces3d = [
                     Meshes.connect((content[p], content[p+1], content[p+2]), Meshes.Triangle) for p = 1:3:length(content)
                 ]
 
                 push!(mesh, "faces" => faces3d)
             elseif i.name == "textureCoords"
-                content = parse_opf_array(i.content)
-                content = Meshes.Point2[
-                    Meshes.Point2(content[[i, i + 1]]) for i in 1:2:length(content)
+                content = parse_opf_array(i.content) ./ 100 * u"m"
+                content = [
+                    Meshes.Point(content[[p, p + 1]]...) for p in 1:2:length(content)
                 ]
                 push!(mesh, "textureCoords" => content)
-            else
-                content = parse_opf_array(i.content)
-                content = Meshes.Point3[
-                    Meshes.Point3(content[[i, i + 1, i + 2]]) for i in 1:3:length(content)
+            elseif i.name == "normals"
+                content = parse_opf_array(i.content) ./ 100 * u"m"
+                content = [
+                    Meshes.Vec(content[[p, p + 1, p + 2]]...) for p in 1:3:length(content)
+                ]
+                push!(mesh, "normals" => content)
+            elseif i.name == "points"
+                content = parse_opf_array(i.content) ./ 100 * u"m"
+                content = [
+                    Meshes.Point(content[[p, p + 1, p + 2]]...) for p in 1:3:length(content)
                 ]
                 push!(mesh, i.name => content)
+            else
+                error("Unknown node element for mesh$(i.Id) in mesh BDD: $(i.name)")
             end
         end
 
@@ -341,7 +349,7 @@ The transformation matrices in `geometry` are 3*4.
 """
 function parse_opf_topology!(node, mtg, features, attr_type, mtg_type, ref_meshes, read_id=true, max_id=Ref(1))
     link = "/" # default, for "topology" and "decomp"
-    b = Vec3(0, 0, 0)
+    b = Vec(0.0u"m", 0.0u"m", 0.0u"m")
     if node.name == "branch"
         link = "+"
     elseif node.name == "follow"
@@ -394,9 +402,9 @@ function parse_opf_topology!(node, mtg, features, attr_type, mtg_type, ref_meshe
                 # See also this for decomposition: https://colab.research.google.com/drive/1ImBB-N6P9zlNMCBH9evHD6tjk0dzvy1_
 
                 #! OK what I could do is use my own transformation function that adds w (=1)
-                #! to the Point3 when transforming it with the 4x4 matrix?
+                #! to the Point when transforming it with the 4x4 matrix?
 
-                transformation = Affine(@view(geom[:mat][1:3, 1:3]), b) → Translate(@view(geom[:mat][1:3, 4])...)
+                transformation = Affine(@view(geom[:mat][1:3, 1:3]), b) → Translate((@view(geom[:mat][1:3, 4]) ./ 100 * u"m")...)
                 # transformation = Translation(geom[:mat][1:3, 4]) ∘ LinearMap(geom[:mat][1:3, 1:3]) # CoordinateTransformations
                 # NB: We read an homogeneous transformation matrix from the OPF, but we work
                 # with cartesian coordinates in PlantGeom by design. So we deconstruct our

src/opf/reference_meshes.jl

@@ -116,7 +116,7 @@ function meshBDD_to_meshes(x)
         mesh_points = pop!(value, "points")
         mesh_faces = pop!(value, "faces")
 
-        points3d = Meshes.Point3[mesh_points[p:(p+2)] for p = 1:3:length(mesh_points)]
+        points3d = Meshes.Point[mesh_points[p:(p+2)] for p = 1:3:length(mesh_points)]
         faces3d = [Meshes.connect((mesh_faces[p:(p+2)]...,), Meshes.Ngon) for p = 1:3:length(mesh_faces)]
 
         push!(mesh, "mesh" => Meshes.SimpleMesh(points3d, faces3d))
@@ -150,17 +150,15 @@ Align all reference meshes along the X axis. Used for visualisation only.
 """
 function align_ref_meshes(meshes::RefMeshes)
     meshes_vec = Meshes.SimpleMesh[]
-    translation_vec = [0.0, 0.0, 0.0]
+    trans = Translate(0.0, 0.0, 0.0)
 
     for i in meshes.meshes
-        translated_vertices = [i + Meshes.Vec(translation_vec...) for i in Meshes.vertices(i.mesh)]
-
-        push!(meshes_vec, Meshes.SimpleMesh(translated_vertices, Meshes.topology(i.mesh)))
-
+        push!(meshes_vec, trans(i.mesh))
         # Maximum X coordinates of the newly translated mesh:
-        xmax_ = maximum([Meshes.coordinates(i)[1] for i in translated_vertices])
+        xmax_ = Meshes.coords(maximum(Meshes.boundingbox(i.mesh))).x
+
         # Update the translation for the next mesh to begin at xmax*1.1 from the last one
-        translation_vec[1] = xmax_ * 1.1
+        trans = Translate(xmax_.val * 1.1, 0.0, 0.0)
     end
 
     return meshes_vec

src/opf/refmesh_to_mesh.jl

@@ -40,7 +40,7 @@ function refmesh_to_mesh(node)
 end
 
 function apply_transformation(transformation, ref_mesh)
-    Meshes.SimpleMesh([transformation(p) for p in ref_mesh.vertices], Meshes.topology(ref_mesh))
+    Meshes.SimpleMesh([transformation(p) for p in Meshes.vertices(ref_mesh)], Meshes.topology(ref_mesh))
 end
 
 function refmesh_to_mesh!(node)