New mesh

Project.toml

@@ -1,10 +1,9 @@
 name = "PlantGeomPrimitives"
 uuid = "7eef3cc5-4580-4ff0-8f6f-933507db6664"
 authors = ["Alejandro Morales Sierra <[email protected]> and contributors"]
-version = "0.0.4"
+version = "0.0.5"
 
 [deps]
-ColorTypes = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
 CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
@@ -14,7 +13,6 @@ Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
-ColorTypes = "0.11.4"
 CoordinateTransformations = "0.6.3"
 FileIO = "1.16.1"
 GeometryBasics = "0.4.7"

src/Mesh/Mesh.jl

@@ -3,24 +3,25 @@
 """
     Mesh
 
-A struct representing a 3D mesh. Every three vertices represents a triangle.
+A struct representing a 3D mesh. Every three vertices represents a triangle. Properties per
+    triangle are stored in a dictionary of arrays.
 
 # Fields
 - `vertices`: A vector containing the vertices of the mesh.
-- `normals`: A vector containing the normals of the mesh.
+- `properties`: A dictionary containing additional properties of the mesh (arrays of properties per triangle).
 
 # Example
 ```jldoctest
 julia> v = [Vec(0.0, 0.0, 0.0), Vec(0.0, 1.0, 0.0), Vec(1.0, 0.0, 0.0)];
 
-julia> n = [Vec(0.0, 0.0, 1.0)];
+julia> p = Dict{Symbol, AbstractVector}(:normal => [Vec(0.0, 0.0, 1.0)]);
 
-julia> m = Mesh(v, n);
+julia> m = Mesh(v, p);
 ```
 """
-struct Mesh{VT}
-    vertices::Vector{VT}
-    normals::Vector{VT}
+struct Mesh{FT}
+    vertices::Vector{Vec{FT}}
+    properties::Dict{Symbol, AbstractVector}
 end
 
 """
@@ -49,7 +50,7 @@ julia> Mesh(Float32);
 ```
 """
 function Mesh(::Type{FT} = Float64) where {FT<:AbstractFloat}
-    Mesh(Vec{FT}[], Vec{FT}[])
+    Mesh(Vec{FT}[], Dict{Symbol, AbstractVector}())
 end
 
 """
@@ -85,7 +86,7 @@ function Mesh(nt::Number, ::Type{FT} = Float64) where {FT<:AbstractFloat}
     sizehint!(v, nv)
     n = Vec{FT}[]
     sizehint!(n, nt)
-    Mesh(v, n)
+    Mesh(v, Dict{Symbol, AbstractVector}(:normal => n))
 end
 
 
@@ -110,12 +111,73 @@ julia> Mesh(verts);
 """
 function Mesh(vertices::Vector{<:Vec})
     VT = eltype(vertices)
-    m = Mesh(vertices, VT[])
+    m = Mesh(vertices, Dict{Symbol, AbstractVector}())
     update_normals!(m)
     return m
 end
 
 
+# Auxilliary function to add properties (from p2 to p1)
+function add_properties!(p1::Dict{Symbol, AbstractVector}, p2::Dict{Symbol, AbstractVector})
+    # Both are empty
+    isempty(p1) && isempty(p2) && (return nothing)
+    # If not, they must have the same properties
+    k1, k2 = (keys(p1), keys(p2))
+    @assert k1 == k2 "Properties of both meshes must be the same"
+    # Add properties (we assume each property is stored in an array-like structure)
+    for k in k1
+        @inbounds append!(p1[k], p2[k])
+    end
+    return p1
+end
+
+"""
+    add_property!(m::Mesh, prop::Symbol, data, nt = ntriangles(m))
+
+Add a property to a mesh. The property is identified by a name (`prop`) and is stored as an
+array of values (`data`), one per triangle. If the property already exists, the new data is
+appended to the existing property, otherwise a new property is created. It is possible to
+pass a single object for `data`, in which case the property will be set to the same value for
+all triangles.
+
+# Arguments
+- `mesh`: The mesh to which the property is to be added.
+- `prop`: The name of the property to be added as a `Symbol`.
+- `data`: The data to be added to the property (an array or a single value).
+- `nt`: The number of triangles to be assumed if `data` is not an array. By default this is the number of triangles in the mesh.
+
+# Returns
+The mesh with updated properties.
+
+# Example
+```jldoctest
+julia> r = Rectangle();
+
+julia> add_property!(r, :absorbed_PAR, [0.0, 0.0]);
+```
+"""
+function add_property!(m::Mesh, prop::Symbol, data, nt = ntriangles(m))
+    # Check if the data is an array and if not convert it to an array with length nt
+    vecdata = data isa AbstractVector ? data : fill(data, nt)
+    # Create new property if the one being added does not exist
+    if !haskey(properties(m), prop)
+        properties(m)[prop] = vecdata
+    # Otherwise append the data to the existing property creating an union type if necessary
+    # when the type of `data` does not inherit from the type of the existing property
+    else
+        etype1 = eltype(properties(m)[prop])
+        etype2 = eltype(vecdata)
+        if etype2 isa etype1
+            append!(properties(m)[prop], vecdata)
+        else
+            etype3 = Union{etype1, etype2}
+            properties(m)[prop] = convert(Vector{etype3}, properties(m)[prop])
+            append!(properties(m)[prop], vecdata)
+        end
+    end
+    return m
+end
+
 """
     Mesh(meshes)
 
@@ -138,57 +200,39 @@ julia> m = Mesh([e,r]);
 """
 function Mesh(meshes::Vector{<:Mesh})
     @assert !isempty(meshes) "At least one mesh must be provided"
-    @inbounds VT = Vec{eltype(first(meshes))}
-    # Positions where each old mesh starts in the new mesh
-    nverts = cumsum(nvertices(m) for m in meshes)
-    nnormals = cumsum(length(normals(m)) for m in meshes)
-    # Allocate elements of the new mesh
-    verts = Vector{VT}(undef, last(nverts))
-    norms = Vector{VT}(undef, last(nnormals))
-    # Fill up the elements of the new mesh
-    @inbounds for i in eachindex(meshes)
-        mesh = meshes[i]
-        # First mesh is simple
-        if i == 1
-            for v = 1:nverts[1]
-                verts[v] = vertices(mesh)[v]
-            end
-            for f = 1:nnormals[1]
-                norms[f] = normals(mesh)[f]
-            end
-        # Other meshes start where previous mesh ended
-        else
-            v0 = nverts[i-1]
-            f0 = nnormals[i-1]
-            for v = v0+1:nverts[i]
-                verts[v] = vertices(mesh)[v-v0]
-            end
-            for f = f0+1:nnormals[i]
-                norms[f] = normals(mesh)[f-f0]
-            end
+    @inbounds verts = copy(vertices(meshes[1]))
+    @inbounds props = properties(meshes[1])
+    if length(meshes) > 1
+        @inbounds for i in 2:length(meshes)
+            append!(verts, vertices(meshes[i]))
+            add_properties!(props, properties(meshes[i]))
         end
     end
-    Mesh(verts, norms)
+    Mesh(verts, props)
 end
 
 # Calculate the normals of a mesh and add them (deals with partially compute normals)
-function update_normals!(m::Mesh)
-    if isempty(m.normals)
-        for i in 1:3:length(m.vertices)
-            @inbounds v1 = m.vertices[i]
-            @inbounds v2 = m.vertices[i+1]
-            @inbounds v3 = m.vertices[i+2]
+function update_normals!(m::Mesh{FT}) where {FT<:AbstractFloat}
+    # 1. Check if there is a property called :normal and if not create it
+    if !haskey(properties(m), :normal)
+        properties(m)[:normal] = Vec{FT}[]
+    end
+    vs = vertices(m)
+    lv = length(vs)
+    # 2. If the property :normal is empty, compute the normals for all vertices
+    if isempty(normals(m))
+        for i in 1:3:lv
+            @inbounds v1, v2, v3 = vs[i], vs[i+1], vs[i+2]
             n = L.normalize(L.cross(v2 .- v1, v3 .- v1))
-            push!(m.normals, n)
+            push!(normals(m), n)
         end
-    elseif length(normals(m)) != div(length(vertices(m)), 3)
-        ln = length(m.normals)
-        for i in 3ln:3:(length(m.vertices) - 3)
-            @inbounds v1 = m.vertices[i + 1]
-            @inbounds v2 = m.vertices[i + 2]
-            @inbounds v3 = m.vertices[i + 3]
+    else
+    # 3. If the property :normal is not empty, compute the normals for the remaining vertices
+        ln = length(normals(m))
+        for i in 3ln:3:(lv - 3)
+            @inbounds v1, v2, v3 = vs[i + 1], vs[i + 2], vs[i + 3]
             n = L.normalize(L.cross(v2 .- v1, v3 .- v1))
-            push!(m.normals, n)
+            push!(normals(m), n)
         end
     end
     return nothing
@@ -206,9 +250,7 @@ Extract the the type used to represent coordinates in a mesh (e.g., `Float64`).
 ```jldoctest
 julia> v = [Vec(0.0, 0.0, 0.0), Vec(0.0, 1.0, 0.0), Vec(1.0, 0.0, 0.0)];
 
-julia> n = [Vec(0.0, 0.0, 1.0)];
-
-julia> m = Mesh(v, n);
+julia> m = Mesh(v);
 
 julia> eltype(m);
 ```
@@ -233,14 +275,12 @@ The number of triangles in the mesh as an integer.
 ```jldoctest
 julia> v = [Vec(0.0, 0.0, 0.0), Vec(0.0, 1.0, 0.0), Vec(1.0, 0.0, 0.0)];
 
-julia> n = [Vec(0.0, 0.0, 1.0)];
-
-julia> m = Mesh(v, n);
+julia> m = Mesh(v);
 
 julia> ntriangles(m);
 ```
 """
-ntriangles(mesh::Mesh) = div(length(mesh.vertices), 3)
+ntriangles(mesh::Mesh) = div(length(vertices(mesh)), 3)
 
 """
     nvertices(mesh)
@@ -257,14 +297,12 @@ The number of vertices in the mesh as an integer.
 ```jldoctest
 julia> v = [Vec(0.0, 0.0, 0.0), Vec(0.0, 1.0, 0.0), Vec(1.0, 0.0, 0.0)];
 
-julia> n = [Vec(0.0, 0.0, 1.0)];
-
-julia> m = Mesh(v, n);
+julia> m = Mesh(v);
 
 julia> nvertices(m);
 ```
 """
-nvertices(mesh::Mesh) = length(mesh.vertices)
+nvertices(mesh::Mesh) = length(vertices(mesh))
 
 # Accessor functions
 """
@@ -282,9 +320,7 @@ A vector containing the vertices of the mesh.
 ```jldoctest
 julia> v = [Vec(0.0, 0.0, 0.0), Vec(0.0, 1.0, 0.0), Vec(1.0, 0.0, 0.0)];
 
-julia> n = [Vec(0.0, 0.0, 1.0)];
-
-julia> m = Mesh(v, n);
+julia> m = Mesh(v);
 
 julia> vertices(m);
 ```
@@ -306,14 +342,36 @@ A vector containing the normals of the mesh.
 ```jldoctest; output=false
 julia> v = [Vec(0.0, 0.0, 0.0), Vec(0.0, 1.0, 0.0), Vec(1.0, 0.0, 0.0)];
 
-julia> n = [Vec(0.0, 0.0, 1.0)];
-
-julia> m = Mesh(v, n);
+julia> m = Mesh(v);
 
 julia> normals(m);
 ```
 """
-normals(mesh::Mesh) = mesh.normals
+normals(mesh::Mesh) = properties(mesh)[:normal]
+
+"""
+    properties(mesh::Mesh)
+
+Retrieve the properties of a mesh. Properties are stored as a dictionary with one entry per
+type of property. Each property is an array of objects, one per triangle. Each property is
+identified by a symbol (e.g.).
+
+# Arguments
+- `mesh`: The mesh from which to retrieve the normals.
+
+# Returns
+A vector containing the normals of the mesh.
+
+# Example
+```jldoctest; output=false
+julia> r = Rectangle();
+
+julia> add_property!(r, :absorbed_PAR, [0.0, 0.0]);
+
+julia> properties(r);
+```
+"""
+properties(mesh::Mesh) = mesh.properties
 
 """
     get_triangle(m::Mesh, i)
@@ -332,9 +390,7 @@ A vector containing the three vertices defining the i-th triangle.
 julia> v = [Vec(0.0, 0.0, 0.0), Vec(0.0, 1.0, 0.0), Vec(1.0, 0.0, 0.0),
             Vec(0.0, 0.0, 0.0), Vec(0.0, 1.0, 0.0), Vec(0.0, 0.0, 1.0)];
 
-julia> n = [Vec(0.0, 0.0, 1.0),Vec(0.0, 0.0, -1.0)];
-
-julia> m = Mesh(v, n);
+julia> m = Mesh(v);
 
 julia> get_triangle(m, 2);
 ```
@@ -413,12 +469,48 @@ function Base.:(==)(m1::Mesh, m2::Mesh)
 end
 
 # Check if two meshes are approximately equal (mostly for testing)
-function Base.isapprox(
-    m1::Mesh,
-    m2::Mesh;
-    atol::Real = 0.0,
-    rtol::Real = atol > 0.0 ? 0.0 : sqrt(eps(1.0)),
-)
+function Base.isapprox(m1::Mesh, m2::Mesh; atol::Real = 0.0,
+                      rtol::Real = atol > 0.0 ? 0.0 : sqrt(eps(1.0)))
     isapprox(vertices(m1), vertices(m2), atol = atol, rtol = rtol) &&
         isapprox(normals(m1), normals(m2), atol = atol, rtol = rtol)
 end
+
+
+"""
+    add!(mesh1, mesh2; kwargs...)
+
+Manually add a mesh to an existing mesh with optional properties captured as keywords. Make
+sure to be consistent with the properties (both meshes should end up with the same lsit of
+properties). For example, if the scene was created with `:colors``, then you should provide
+`:colors`` for the new mesh as well.
+
+# Arguments
+- `mesh1`: The current mesh we want to extend.
+- `mesh1`: A new mesh we want to add.
+- `kwargs`: Properties to be set per triangle in the new mesh.
+
+# Example
+```jldoctest
+julia> t1 = Triangle(length = 1.0, width = 1.0);
+
+julia> using ColorTypes: RGB
+
+julia> add_property!(t1, :colors, rand(RGB));
+
+julia> t2 = Rectangle(length = 5.0, width = 0.5);
+
+julia> add!(t1, t2, colors = rand(RGB));
+```
+"""
+function add!(mesh1, mesh2; kwargs...)
+    # Make sure the mesh contains normals
+    update_normals!(mesh2)
+    # Add the vertices and normals
+    append!(vertices(mesh1), vertices(mesh2))
+    add_property!(mesh1, :normal, normals(mesh2))
+    # Set optional properties per triangle
+    for (k, v) in kwargs
+        add_property!(mesh1, k, v, ntriangles(mesh2))
+    end
+    return mesh1
+end