Add implementation of ellipsoid

Author: AleMorales

NEWS.md

@@ -2,6 +2,12 @@
 
 We started keeping track of changes in the `NEWS.md` file after version 0.1.0.
 
+# PlantGeomPrimitives 1.0.1 (2026-06-16)
+
+Implement `Ellipsoid` as a primitive discretize into rectangular patches defined by fixed longitude and 
+latitude increments. The user specifies the number of longitude and latitude increments rather than the
+total number of triangles.
+
 # PlantGeomPrimitives 1.0.1 (2026-06-11)
 
 The rotation functions were still expecting angles in radians rather than degrees. The

Project.toml

@@ -1,7 +1,7 @@
 name = "PlantGeomPrimitives"
 uuid = "7eef3cc5-4580-4ff0-8f6f-933507db6664"
 authors = ["Alejandro Morales Sierra <alejandro.moralessierra@wur.nl> and contributors"]
-version = "1.0.1"
+version = "1.0.2"
 
 [deps]
 ColorTypes = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"

src/Primitives/Ellipsoid.jl

@@ -1,61 +1,106 @@
 ### This file contains public API ###
 
-function equal_solid_angles(n)
-    # Distribution sectors along zenith angles
-    Δθ = π / 2 / n
-    uθₗ = 0.0:Δθ:π/2-Δθ
-    uθᵤ = Δθ:Δθ:π/2
-
-    # Calculate number of azimuth sectors and ΔΦ per zenith ring
-    fac = cos.(uθₗ) - cos.(uθᵤ)
-    n2 = n * n
-    ns = round.(fac ./ sum(fac) .* n2)
-    sum(ns) != n2 && (ns[end] = ns[end] + n2 - sum(ns))
-    ΔΦs = 2π ./ ns
-
-    # Generate coordinates of all sectors
-    c = 1
-    θₗ = Vector{Float32}(undef, n2)
-    θᵤ = Vector{Float32}(undef, n2)
-    Φₗ = Vector{Float32}(undef, n2)
-    Φᵤ = Vector{Float32}(undef, n2)
-    for i = 1:n
-        ΔΦ = ΔΦs[i]
-        for j = 1:ns[i]
-            θₗ[c] = Δθ * (i - 1)
-            θᵤ[c] = Δθ * i
-            Φₗ[c] = ΔΦ * (j - 1)
-            Φᵤ[c] = ΔΦ * j
-            c += 1
+function vertex_sphere(θ, φ, trans)
+    FT = eltype(trans.linear)
+    sθ, cθ = sincos(FT(θ))
+    sφ, cφ = sincos(FT(φ))
+    trans(Vec{FT}(sθ * cφ, sθ * sφ, 1 + cθ))
+end
+
+struct EllipsoidVertices{FT,TT}
+    n::Int
+    dθ::FT
+    dφ::FT
+    trans::TT
+end
+function EllipsoidVertices(n, trans)
+    FT = eltype(trans.linear)
+    EllipsoidVertices(n, FT(π / n), FT(2π / n), trans)
+end
+
+function Base.iterate(c::EllipsoidVertices{FT,TT}, i::Int = 1)::Union{Nothing,Tuple{Vec{FT},Int64}} where {FT,TT}
+    n  = c.n
+    i > 6n * (n - 1) && return nothing
+    ntop = 3n
+    nb   = 6n  # vertices per middle band
+    if i <= ntop
+        # Top cap: n fan triangles from north pole to first latitude ring
+        j = div(i - 1, 3)
+        v = mod(i - 1, 3) + 1
+        v == 1 && (return c.trans(Vec{FT}(0, 0, 2)), i + 1)
+        v == 2 && (return vertex_sphere(c.dθ,       j       * c.dφ, c.trans), i + 1)
+        v == 3 && (return vertex_sphere(c.dθ,      (j + 1)  * c.dφ, c.trans), i + 1)
+    elseif i <= ntop + (n - 2) * nb
+        # Middle bands: (n-2) latitude bands, each with n quads split into 2 triangles
+        im    = i - ntop - 1
+        iband = div(im, nb)
+        iw    = mod(im, nb)
+        θu    = (iband + 1) * c.dθ   # upper ring (smaller θ)
+        θl    = (iband + 2) * c.dθ   # lower ring (larger θ)
+        if iw < 3n
+            # First triangle of each quad: (top-left, bottom-left, bottom-right)
+            j = div(iw, 3)
+            v = mod(iw, 3) + 1
+            v == 1 && (return vertex_sphere(θu,  j      * c.dφ, c.trans), i + 1)
+            v == 2 && (return vertex_sphere(θl,  j      * c.dφ, c.trans), i + 1)
+            v == 3 && (return vertex_sphere(θl, (j + 1) * c.dφ, c.trans), i + 1)
+        else
+            # Second triangle of each quad: (top-left, bottom-right, top-right)
+            iw2 = iw - 3n
+            j   = div(iw2, 3)
+            v   = mod(iw2, 3) + 1
+            v == 1 && (return vertex_sphere(θu,  j      * c.dφ, c.trans), i + 1)
+            v == 2 && (return vertex_sphere(θl, (j + 1) * c.dφ, c.trans), i + 1)
+            v == 3 && (return vertex_sphere(θu, (j + 1) * c.dφ, c.trans), i + 1)
         end
+    else
+        # Bottom cap: n fan triangles from last latitude ring to south pole (reversed winding)
+        ib    = i - ntop - (n - 2) * nb - 1
+        j     = div(ib, 3)
+        v     = mod(ib, 3) + 1
+        θlast = (n - 1) * c.dθ
+        v == 1 && (return c.trans(Vec{FT}(0, 0, 0)), i + 1)
+        v == 2 && (return vertex_sphere(θlast, (j + 1) * c.dφ, c.trans), i + 1)
+        v == 3 && (return vertex_sphere(θlast,  j      * c.dφ, c.trans), i + 1)
     end
-
-    (θₗ, θᵤ, Φₗ, Φᵤ, ns)
 end
 
-# Every combination of theta and phi angles yields a vertex
-# Then the vertex are mapped to the other hemisphere
-# This means that we cannot have vertices at theta = 0
+Base.length(c::EllipsoidVertices) = 6c.n * (c.n - 1)
+Base.eltype(::Type{EllipsoidVertices{FT,TT}}) where {FT,TT} = Vec{FT}
 
-# Step 1: Turn the function above into an interator that generates
-# vertices in both hemispheres (following a reasonable ordering)
 
-# Step 2: Create the rules to connect these vertices into triangles
+"""
+    Ellipsoid(;length = 1.0, width = 1.0, height = 1.0, n = 20)
 
-# Step 3: Create the normals to each face (average theta and phi??)
+Create a solid ellipsoid with dimensions given by `length`, `width` and `height`,
+with the base on the z = 0 plane, discretized using a UV latitude-longitude mesh
+with `n` latitude and azimuth subdivisions.
 
+# Arguments
+- `length = 1.0`: The length of the ellipsoid along the z axis.
+- `width = 1.0`: The width of the ellipsoid along the y axis.
+- `height = 1.0`: The height of the ellipsoid along the x axis.
+- `n = 20`: Resolution (number of latitude and azimuth subdivisions). Must be >= 3.
+  The mesh has `2n(n-1)` triangles.
 
-# Scaled ellipsoid
-function Ellipsoid(l::Number, w::Number, h::Number, n::Number)
-    @error "Ellipsoid not implemented yet"
+# Examples
+```jldoctest
+julia> Ellipsoid(;length = 1.0, width = 1.0, height = 1.0, n = 20);
+```
+"""
+function Ellipsoid(; length::FT = 1.0, width::FT = 1.0, height::FT = 1.0, n::Int = 20) where {FT}
+    trans = CT.LinearMap(CT.SDiagonal(height / FT(2), width / FT(2), length / FT(2)))
+    Ellipsoid(trans; n = n)
 end
 
-# Create a ellipsoid from affine transformation
-function Ellipsoid(trans::CT.AbstractAffineMap, nt::Number)
-    @error "Ellipsoid not implemented yet"
+# Create an ellipsoid from affine transformation
+function Ellipsoid(trans::CT.AbstractAffineMap; n::Int = 20)
+    @assert n >= 3 "n must be at least 3"
+    Primitive(trans, x -> EllipsoidVertices(n, x))
 end
 
-# Create a ellipsoid from affine transformation and add it in-place to existing mesh
-function Ellipsoid(m::Mesh, trans::CT.AbstractAffineMap, nt::Int)
-    @error "Ellipsoid not implemented yet"
+# Create an ellipsoid from affine transformation and add it in-place to existing mesh
+function Ellipsoid!(m::Mesh, trans::CT.AbstractAffineMap; n::Int = 20)
+    @assert n >= 3 "n must be at least 3"
+    Primitive!(m, trans, x -> EllipsoidVertices(n, x))
 end

test/runtests.jl

@@ -58,6 +58,9 @@ end
 @testset "solid_cone" begin
     include("test_solid_cone.jl")
 end
+@testset "ellipsoid" begin
+    include("test_ellipsoid.jl")
+end
 @testset "transformations" begin
     include("test_transformations.jl")
 end

test/test_ellipsoid.jl

@@ -3,7 +3,70 @@ using Test
 import CoordinateTransformations: SDiagonal, LinearMap
 
 let
+    n = 20
 
+    # Standard Ellipsoid primitive
+    e = G.Ellipsoid(length = 2.0, width = 1.0, height = 1.0, n = n)
+    @test e isa G.Mesh
+    @test G.nvertices(e) == 6n * (n - 1)
+    @test G.ntriangles(e) == div(G.nvertices(e), 3)
+    @test length(G.normals(e)) == G.ntriangles(e)
 
+    # Area test: sphere of diameter 2 (length=width=height=2) -> area = 4π
+    e_sphere = G.Ellipsoid(length = 2.0, width = 2.0, height = 2.0, n = n)
+    @test abs(G.area(e_sphere) / (4π) - 1.0) < 0.02
 
+    # Outward normals: each face normal must point away from the ellipsoid center
+    center = G.Vec(0.0, 0.0, 1.0)  # canonical (0,0,1) mapped by SDiagonal(1,1,1)
+    for i in 1:G.ntriangles(e_sphere)
+        tri      = G.get_triangle(e_sphere, i)
+        centroid = (tri[1] .+ tri[2] .+ tri[3]) ./ 3
+        @test sum(G.normals(e_sphere)[i] .* (centroid .- center)) > 0
+    end
+
+    # Float32 variant
+    e32 = G.Ellipsoid(length = 2.0f0, width = 1.0f0, height = 1.0f0, n = n)
+    @test e32 isa G.Mesh
+    @test G.nvertices(e32) == 6n * (n - 1)
+    @test G.ntriangles(e32) == div(G.nvertices(e32), 3)
+    @test length(G.normals(e32)) == G.ntriangles(e32)
+
+    # Merging two meshes
+    e2 = G.Ellipsoid(length = 3.0, width = 0.5, height = 0.5, n = n)
+    function foo()
+        e  = G.Ellipsoid(length = 2.0, width = 1.0, height = 1.0, n = n)
+        e2 = G.Ellipsoid(length = 3.0, width = 0.5, height = 0.5, n = n)
+        G.Mesh([e, e2])
+    end
+    m = foo()
+    @test G.nvertices(m)  == G.nvertices(e)  + G.nvertices(e2)
+    @test G.ntriangles(m) == G.ntriangles(e) + G.ntriangles(e2)
+    @test abs(G.area(m) - (G.area(e) + G.area(e2))) < 1e-10
+
+    # Affine map equivalence: keyword constructor and explicit LinearMap must agree
+    scale = LinearMap(SDiagonal(0.5 / 2, 0.5 / 2, 0.5 / 2))
+    e3 = G.Ellipsoid(scale; n = n)
+    e4 = G.Ellipsoid(length = 0.5, width = 0.5, height = 0.5, n = n)
+    @test G.normals(e3)   == G.normals(e4)
+    @test G.vertices(e3)  == G.vertices(e4)
+
+    # In-place Ellipsoid!
+    function foo2()
+        m2     = G.Ellipsoid(length = 2.0, width = 1.0, height = 1.0, n = n)
+        scale2 = LinearMap(SDiagonal(0.25 / 2, 0.5 / 2, 1.5 / 2))
+        G.Ellipsoid!(m2, scale2; n = n)
+        m2
+    end
+    m2 = foo2()
+    e5 = G.Ellipsoid(length = 1.5, width = 0.5, height = 0.25, n = n)
+    @test G.nvertices(m2)  == G.nvertices(e) + G.nvertices(e5)
+    @test G.ntriangles(m2) == G.ntriangles(e) + G.ntriangles(e5)
 end
+
+
+#= import GLMakie
+import PlantViz as PV
+e = G.Ellipsoid(length = 1.0, width = 1.0, height = 1.0, n = n)
+G.add_property!(e, :colors, RGB(1.0, 0.0, 0.0))
+PV.render(e, normals = true)
+=#