Procedural Plant Modeling ,Computer-Aided Design ,Volumetric Mesh Generation ,Finite Element Analysis ,Dynamic Plant Modeling

1 Introduction↩︎

Accurate geometric modeling of plant architecture is a fundamental requirement in computational agricultural applications, such as digital twins [1]–[4], functional-structural plant models [5]–[8], precision agriculture [9], [10], resource optimization [11], [12], genotype-to-phenotype modeling [13]–[15], bio-design [16]–[18], and crop breeding [19]–[22]. High-fidelity geometric models generated from advanced sensor datasets provide quantitative analysis of plant architecture traits under different genetic, environmental, and management conditions, improving agricultural productivity and resilience [8], [12], [23], [24]. By integrating real-time sensor data and environmental parameters, plant models can provide a realistic virtual representation of plants, enabling non-invasive, dynamic, and bi-directional decision-making. In precision agriculture, they allow for accurate morphological evaluations, canopy architecture design, efficient plant breeding, and resource usage [20], [25]. In simulation-driven crop design, high-fidelity crop models can be used in inverse design applications to generate idealized ideotypes for high productivity and improved resilience [26]–[30]. The goal is not only to generate plant models for visualization but also to use them directly for physics-based simulations, such as nutrient and water uptake, structural mechanics [31]–[34], light-use efficiency, and ecophysiological and ecosystem modeling [35]–[41].

The advent of advanced sensor technologies for three-dimensional (3D) plant reconstruction, such as light detection and ranging (LiDAR), terrestrial laser scanning (TLS), multiview stereo (MVS), structured light scanning, neural radiance fields (NeRF), and Gaussian splatting, has enabled non-invasive, high-throughput plant phenotyping and monitoring [37], [42]–[47]. Whole plant architecture is typically represented as depth maps, point clouds, or surface meshes. These reconstructions provide high-resolution quantitative and qualitative visual representations of plant architecture, but they are not directly suitable for physics-based simulations [19], [48]–[52]. This identifies a critical research gap, where we have access to large datasets of high-fidelity 3D plant reconstructions acquired from sensors. However, converting them into analysis-ready volumetric meshes remains computationally intensive, requires manual effort, and is not scalable for complex geometries and challenging applications.

Procedural modeling approaches encode plant geometry in a rule-based manner, mathematically embedding the latent plant branching architecture and topology. Grammar-based approaches, such as L-systems, encode branching topology and plant growth using string-based rules [53]–[60]. Functional-structural plant models (FSPMs) utilize algorithmic rules to construct plant geometry and employ physiological processes to guide the development of plant morphology [5], [7], [38], [61]–[65]. Non-uniform rational B-splines (NURBS) have been used in procedural modeling frameworks to accurately generate high-fidelity plant representations [19], [50], [66]–[68]. These approaches allow for biologically realistic plant modeling, but they are limited to implicit or surface representations, which require considerable computational effort to generate analysis-suitable volumetric meshes.

Generating analysis-suitable volumetric meshes for branching structures is an especially challenging and non-trivial problem, particularly near junctions, where multiple branch segments connect at different angles. A robust volumetric modeling framework should capture the correct plant branching topology for all possible configurations while generating a regular and conforming mesh that is free of intersections, gaps, and distortions [69], [70]. Conventional volumetric meshing and NURBS-based modeling workflows for generating branching geometries can create highly accurate and high-quality models [71]–[75]. However, these approaches are static; modifying the branching topology, such as adding or deleting branches or changing the branch radius, requires considerable effort to remesh or even rebuild the entire geometry from scratch. This highlights a key limitation: there are currently no procedural volumetric meshing frameworks for automatically generating branching geometries while remaining compatible with changing plant topology.

We present an automated volumetric procedural modeling framework for generating analysis-suitable hexahedral (hex) meshes of plant branching structures from 3D point clouds and their corresponding plant skeletons [76]–[82]. The input skeleton is first converted into a rooted directional graph representation that captures the plant topology. B-spline curves are then fitted to each branch, rendering a smooth centerline geometry. Rotation-minimizing frames along the centerline curves are used to construct cylindrical B-spline volumes for individual branches. At bifurcation and trifurcation junctions, incident branch B-spline volumes are joined via control-point blending operations, enforcing seamless transitions. The resulting volumetric parameterization of the complete plant is then evaluated to generate a conforming and smooth hex mesh suitable for finite element analysis (FEA).

The proposed framework utilizes a compact set of rule-based generation steps instead of manual mesh construction. It generates plant geometry through a sequence of steps, progressing from graph and curve representations to volumetric meshes. The framework not only enables automatic volumetric mesh generation from plant skeleton data but also allows local updates without reconstructing the whole plant geometry as the plant topology changes. The main contributions of our work are:

1. A procedural volumetric modeling framework that automatically converts 3D plant skeleton data obtained from sensors into hex meshes suitable for FEA.

2. A branch-wise B-spline volumetric parameterization strategy that first represents each branch centerline as a B-spline curve. Rotation-minimizing frames are constructed along the centerline curves, and tensor-product cylindrical B-spline volumes are generated for each branch.

3. A junction-blending procedure is applied to bifurcations and trifurcations, joining incident branch volumes with seamless transitions and producing conforming volumetric meshes for the entire plant geometry.

4. Demonstration of analysis-suitable volumetric mesh generation for diverse plant species of increasing geometric complexity, which includes mung bean, tomato, and walnut trees.

5. Robust generation of plant hex meshes with spatially varying branch radii from unstructured 3D point clouds.

6. Demonstration of dynamic mesh generation for plant growth modeling, where newly added branches are locally parameterized and blended without reconstructing the whole plant geometry.

2 Related work↩︎

In recent years, there has been significant advancement in the field of 3D plant reconstruction technologies, enabling the capture of detailed structural and morphological plant traits for accurate phenotyping. These include active depth-based imaging technologies such as LiDAR [83], TLS [84]–[87], structured-light sensors [88], [89], and time-of-flight cameras [90], which use light or infrared signals to generate depth maps. Passive 3D reconstruction methods, such as stereo vision, MVS [91]–[94], and structure-from-motion(SfM) photogrammetry [95], offer ways to reconstruct depth maps without emitting signals. However, these methods generate point clouds that suffer from noise, occlusions, and missing data [23]. More recently, deep learning-based approaches such as NeRF [96], [97] and 3D Gaussian splatting [98] have emerged as phenotyping tools that generate high-resolution point clouds [49], [99]–[102]. These approaches generate continuous volumetric scene representations, allowing them to handle thin structures and occlusions more easily [42], [103]. After 3D point cloud acquisition, the phenotyping process proceeds with image analysis steps, such as point cloud registration [104], surface reconstruction [19], [67], [105], [106], semantic segmentation [107], and skeletonization [108], [109], to extract key morphological traits [46], [89], yielding segmented point clouds, voxels, surfaces, or skeletons that describe plant geometry. Current approaches focus primarily on surface reconstruction for plant phenotyping, using triangulations [110], [111], B-splines, or NURBS [112], [113] to represent plant architecture at the organ, plant, or canopy scales [50], [106]. As a result, analysis-suitable volumetric modeling remains largely absent from plant phenotyping.

Procedural plant modeling frameworks utilize algorithmic rules rather than manual geometry generation to model 3D plant architecture. Using a small set of geometric parameters and a rule-based generation approach, geometric models can be generated more efficiently for realistic plant growth and phenotyping applications [114], [115]. Lindenmayer Systems (L-systems) are rule-based procedural modeling approaches used to generate realistic 3D plant structures through a recursive algorithmic generation of L-strings that encode plant topology, geometry, and growth [116]. Deterministic, stochastic, and parametric L-systems have been shown to model realistic 1D plant branching architectures and growth processes for plant visualization and botany [54], [117], [118]. Subsequent research extended these approaches to include interactive and sketch-based tree modeling and procedural reconstruction from images or point clouds [57], [119]–[121]. These approaches showcase the flexibility of the rule-based frameworks for both forward and inverse data-driven reconstruction. FSPM approaches utilize physiological processes along with plant geometry to model plant growth and function [7], [63], [122]. The geometric modeling components, such as plant organs and internodes, are combined with mechanistic modules for photosynthesis, carbon allocation, and resource competition in a unified simulation framework [8], [38], [64], [87], [123]. Both the L-system and FSPMs provide geometric outputs for visualization, trait analysis, and growth simulation; however, they cannot be directly used for physics-based simulations [63], [122].

Recently, data-driven approaches have utilized machine learning to encode key morphological parameters, such as plant topology and geometry, into compact representations [102], [124]–[126]. By training on large imaging and point cloud datasets, supervised and deep learning approaches, such as convolutional neural networks (CNNs) and multi-layer perceptrons (MLPs), recover parameters for morphological traits and plant growth, enabling automated trait estimation and complementing rule-based or mechanistic models [126]. They also provide low-dimensional embeddings of complex plant topology and geometry for inverse modeling and generative modeling of diverse plant geometries under varying environmental and genetic conditions. Computer-aided design (CAD) modeling allows for a smooth parametric representation of complex geometry. Spline-based modeling allows for precise control, enabling both local and global high-fidelity shape representation [113]. B-spline curves or generalized cylinders have been used to represent plant branches, crowns, roots, and plant organs [127]–[129]. A generalized cylinder-based approach for modeling plant branches with implicit ramiform junctions was introduced in [130]. In [131], plant branches were represented as 3D B-spline curves coupled with a 1D B-spline radius; however, junction blending was identified as an open problem. Generalized B-spline crown envelopes have been used to model the complex polymorphism of Chinese fir to support natural pruning simulations [132]. B-splines have also been coupled with L-systems for the visual modeling of rice root growth [133], [134]. Recent methods utilize B-spline centerline approaches along with parallel transport frames for plant geometric modeling. Procedural frameworks for NURBS-based surface fitting to point cloud data have been developed, with clear gaps remaining in analysis-suitable volumetric modeling [19], [50], [65], [135]–[137].

Volumetric hex meshing and spline-based volumetric parameterization (V-rep) methods have been developed for FEA and isogeometric analysis (IGA) for branching structures [73], [74], [138]–[149]. Trivariate tensor-product B-spline volumes were fitted to genus-0 surfaces using harmonic parameterization in the cylindrical coordinate system in [150]. In [151], trimmed tri-variate B-spline volumes using geometric primitives were glued along shared boundaries with \(C^{0}\)-continuity. Patient-specific NURBS volumetric modeling has been carried out for vascular structures using skeleton-based sweeps with template-based junction continuity enforcement [72]. Polycube and block-decomposition approaches have been used to generate high-quality, analysis-suitable volumetric models for complex geometries [152]–[159]. These robust, analysis-suitable modeling approaches, however, are not tailored for plant phenotyping workflows and typically assume static geometric construction, with limited applicability for modeling dynamic plant branching topologies.

3 Methods↩︎

We present a procedural framework that generates hex meshes for FEA from input plant skeleton point clouds. These point clouds approximate the medial axis of the plant, capturing the branching topology and centerline geometry. We denote the 3D point cloud by \(\mathbf{X}=\{\mathbf{x}_l\}_{l=1}^{N}\), with \(\mathbf{x}_l\in\mathbb{R}^3\), where \(N\) denotes the total number of points. The input point cloud is first converted into a skeleton graph that encodes the branching topology and connectivity. The skeleton graph is then transformed into an analysis-suitable hex mesh through five steps: (1) extraction of a branch graph from the skeleton graph by collapsing maximal chains of degree-two vertices; (2) fitting a B-spline centerline curve to the coordinate sequence along each edge of the branch graph; (3) construction of rotation-minimizing reference frames along each branch centerline curve by parallel transport; (4) construction of a hollow cylindrical tensor-product B-spline volume for each branch; and (5) blending of control points at each junction node. The B-spline-based volumetric parameterization provides a smooth, locally controllable representation of the branch geometry. The compact support of B-spline basis functions confines junction blending to a small number of control point layers in the volume. The output of the framework is a conforming hex mesh obtained by evaluating and assembling the blended B-spline branch volumes. The detailed procedure for each of these steps is described in the remainder of this section.

3.1 Skeleton and branch graph representation↩︎

Let \(\mathcal{G}_S = (\mathcal{V}_S, \mathcal{E}_S)\) denote the skeleton graph, with nodes \(\nu_l \in \mathcal{V}_S\) and edges \(e_{m} \in \mathcal{E}_S\). We assume that the input skeleton graph has been cleaned to form a connected graph that captures the correct branching topology; skeleton extraction and cleanup from the raw point cloud are outside the scope of this work. Each node \(\nu_l\) is associated with an embedded coordinate \(\mathbf{x}_l\) and a non-negative scalar \(R_{l} \ge 0\) representing the outer radius of the cross-section at that coordinate. From \(\mathcal{G}_S\), we then construct a branch graph \(\mathcal{G}_B = (\mathcal{V}_B, \mathcal{E}_B)\) that is treated as a rooted tree, with all edges oriented from the root toward the leaves. The node set \(\mathcal{V}_B \subseteq \mathcal{V}_S\) is defined by retaining all nodes in \(\mathcal{V}_S\) whose degree is not equal to two: \[\mathcal{V}_B = \{\nu_l \in \mathcal{V}_S \mid \deg_{\mathcal{G}_S}(\nu_l) \neq 2\}.\] Here, the degree of the node, \(\deg_{\mathcal{G}_S}(\nu_l)\), is the number of edges in \(\mathcal{G}_S\) incident to \(\nu_l\). Nodes with \(\deg_{\mathcal{G}_S}(\nu_l) = 1\) are terminal nodes. Among the terminal nodes, we select the root node, denoted by \(\nu_{\mathrm{root}}\), with \(\mathcal{V}^{\mathrm{root}}_B=\{\nu_{\mathrm{root}}\}\) associated with either the basal skeleton vertex of the point cloud, such as the vertex corresponding to the lowest \(z\)-coordinate, or a user-defined vertex. The remaining terminal nodes form the leaf set \(\mathcal{V}^{\mathrm{leaf}}_B\). Nodes with \(\deg_{\mathcal{G}_S}(\nu_l) \geq 3\) are the junction nodes. Every junction node in the junction set \(\mathcal{V}^{\mathrm{junc}}_B\) has one incoming edge and at least two outgoing edges in \(\mathcal{G}_B\). The node set \(\mathcal{V}_B\) can thus be partitioned into the root, junction, and leaf sets as \[\mathcal{V}_B = \mathcal{V}^{\mathrm{root}}_B \;\dot{\cup}\; \mathcal{V}^{\mathrm{junc}}_B \;\dot{\cup}\; \mathcal{V}^{\mathrm{leaf}}_B. \label{eq:node95partition}\tag{1}\] Each directed edge from a proximal node \(\nu_\beta\) to a distal node \(\nu_{\alpha}\) is denoted by \(b_{\alpha} \in \mathcal{E}_B\). We identify each edge \(b_{\alpha}\) by its distal node \(\nu_{\alpha}\). Each edge \(b_{\alpha}\) stores an ordered coordinate sequence \(\mathbf{Q}^{\alpha} = (\mathbf{q}^{\alpha}_1, \dots, \mathbf{q}^{\alpha}_{M^{\alpha}})\), corresponding to the unique maximal path in \(\mathcal{G}_S\) from \(\nu_\beta\) to \(\nu_{\alpha}\), whose interior nodes all have degrees equal to two. The endpoints satisfy \(\mathbf{q}^{\alpha}_1=\mathbf{x}_{\beta}\) and \(\mathbf{q}^{\alpha}_{M^{\alpha}}=\mathbf{x}_{\alpha}\), where \(M^{\alpha}\) denotes the number of skeleton points along each edge \(b_{\alpha}\).

3.2 B-spline centerline curve construction↩︎

We represent the centerline geometry of each edge \(b_{\alpha} \in \mathcal{E}_B\) by a smooth centerline B-spline curve, \(\mathbf{C}^{\alpha}:[0,1]\rightarrow\mathbb{R}^3\), defined as: \[\mathbf{C}^\alpha(u) \;=\; \sum_{k=1}^{n^\alpha} N_{k,p}(u)\, \mathbf{P}^\alpha_k, \qquad u\in[0,1], \label{eq:centerline}\tag{2}\] where \(n^\alpha\) denotes the number of control points, \(N_{k,p}(u)\) denotes the \(k^{th}\) univariate basis function of degree \(p\) defined on the clamped knot vector \(\mathcal{U}^{\alpha} := \{u_1, u_2, \dots, u_{n^{\alpha} + p + 1}\}\), and \(\mathbf{P}^\alpha_k \in \mathbb{R}^3\) is the corresponding B-spline control point [160]. We present two strategies to construct the centerline curves for the entire branch graph. The global construction strategy (Section 3.2.1) extracts the complete branch graph and fits an individual B-spline curve to each edge \(b_{\alpha}\). In the hierarchical construction strategy (Section 3.2.2), the branch graph is initialized with a primary path (main stem) in the skeleton graph from the root node to a terminal leaf node. A B-spline curve is first fitted to this primary path to represent the centerline geometry. The complete branch graph is then grown recursively by attaching child branches one at a time. At each new junction, a knot insertion operation is performed on the B-spline curve to split the parent branch, create space for attaching the child branch, and enforce junction blending. In both strategies, the B-spline curve \(\mathbf{C}^\alpha\) is fitted to a target coordinate sequence \(\mathbf{D}^{\alpha}=(\mathbf{d}^{\alpha}_1, \dots, \mathbf{d}^{\alpha}_{M^{\alpha}_{D}})\) of edge \(b_\alpha\), consisting of \(M_D^{\alpha}\) ordered coordinates. In the global construction strategy, we set \(\mathbf{D}^{\alpha}=\tilde{\mathbf{Q}}^{\alpha}\), where \(\tilde{\mathbf{Q}}^{\alpha}\) denotes the coordinate sequence obtained by pruning \(\mathbf{Q}^{\alpha}\) to remove coordinates near each junction node (Section 3.2.1); in the hierarchical construction strategy, we instead set \(\mathbf{D}^{\alpha}=\mathbf{Q}^{\alpha}\). For each coordinate \(\mathbf{d}^{\alpha}_\kappa\), \(\kappa = 1, \dots, M_D^{\alpha}\), we first compute the corresponding B-spline parameter \(\bar{u}_\kappa\) by chord-length parameterization [160]: \[\bar{u}_\kappa \;=\; \frac{\sum_{j=2}^{\kappa} \|\mathbf{d}_j^{\alpha} - \mathbf{d}_{j-1}^{\alpha}\|}{\sum_{j=2}^{M^{\alpha}_D} \|\mathbf{d}_j^{\alpha} - \mathbf{d}_{j-1}^{\alpha}\|}, \qquad \bar{u}_1 = 0, \quad \bar{u}_{M^{\alpha}_D} = 1. \label{eq:chord95length}\tag{3}\] The fitting loss is defined as the mean squared distance between the B-spline curve \(\mathbf{C}^{\alpha}\) and the target coordinate sequence \(\mathbf{D}^{\alpha}\): \[\mathcal{L}^{\alpha}\big(\mathbf{C}^{\alpha}\big) \;=\; \frac{1}{M^{\alpha}_D} \sum_{\kappa=1}^{M^{\alpha}_D} \big\|\, \mathbf{C}^{\alpha}(\bar{u}_\kappa) - \mathbf{d}_\kappa^{\alpha} \,\big\|^2. \label{eq:fit95loss}\tag{4}\] In both construction strategies, the B-spline control point positions, \(\mathbf{P}^\alpha\), are optimized to minimize the loss function \(\mathcal{L}^{\alpha}\) using Adam optimization [161]. Fig. 1 summarizes the complete pipeline from skeleton graph construction to B-spline centerline curve fitting.

3.2.1 Global construction strategy↩︎

In the global construction strategy, the entire branch graph \(\mathcal{G}_B\) is extracted from \(\mathcal{G}_S\). Before B-spline curve fitting, the ordered coordinate sequence \(\mathbf{Q}^{\alpha}\) is pruned near each junction node to leave room for subsequent junction blending between adjacent B-spline volumes (Section 3.4). At each junction node \(\nu_\delta\in\mathcal{V}_B^{\mathrm{junc}}\), we define a pruning sphere of radius \(R^{\mathrm{prune}}_{\delta}\) centered at its embedded coordinate \(\mathbf{x}_\delta\). \(R^{\mathrm{prune}}_{\delta}\) is set as the largest radius of the cross-section of all the branches incident at \(\nu_\delta\). Any coordinate of \(\mathbf{Q}^{\alpha}\) that lies inside this sphere is removed, and the surviving ordered coordinate sequence forms the target coordinate set \(\tilde{\mathbf{Q}}^{\alpha}\) used for B-spline curve fitting. The detailed steps for the pruning procedure are outlined in Algorithm 2.

For each edge \(b_\alpha\), the number of control points for each centerline B-spline curve is assigned based on its branch length, with longer branches assigned more control points. The parameters of the B-spline curve are calculated for the target sequence \(\tilde{\mathbf{Q}}^{\alpha}\) using Eq. 3 . Finally, for each edge \(b_\alpha\), the B-spline curve \(\mathbf{C}^{\alpha}\) is fitted to \(\tilde{\mathbf{Q}}^{\alpha}\) by minimizing the loss \(\mathcal{L}^{\alpha}\) with respect to \(\mathbf{P}^{\alpha}\) (Eq. 4 ). In this construction strategy, an independent B-spline curve is fitted to each branch, so adjacent branch curves meeting at the junction do not share a common tangent direction. Tangential continuity is therefore not guaranteed and is enforced separately during junction blending between the branch volumes (Section 3.4). Algorithm 3 outlines the steps carried out in the global construction strategy.

3.2.2 Hierarchical construction strategy↩︎

In the hierarchical construction strategy, we construct the complete branch graph recursively, starting from a primary path and progressively attaching child branches according to the proximal-to-distal branching order. A primary path, referred to as the main stem, is identified as the path from the root node to a selected terminal leaf node in \(\mathcal{G}_S\). The branch graph \(\mathcal{G}_B\) for the main stem is initialized with the root node \(\nu_{\mathrm{root}}\), the terminal leaf node \(\nu_{\alpha}\), and a single directed edge \(b_{\alpha} \in \mathcal{E}_B\) connecting them. By default, the terminal leaf node \(\nu_{\alpha}\) is selected as the leaf node at the end of the longest root-to-leaf path; although a user-defined leaf node can also be selected. Let \(\mathbf{Q}^{\alpha}\) denote the ordered coordinate sequence along the edge \(b_{\alpha}\). The coordinates in \(\mathbf{Q}^{\alpha}\) are used as target coordinates for B-spline curve fitting (Eqs. 3 4 ). Unlike the global construction strategy, \(n^{\alpha}\) for each edge is not set according to branch length-based binning, as not all the branch edges are available at initialization. Instead, \(n^{\alpha}\) is set equal to the number of coordinates in \(\mathbf{Q}^{\alpha}\).

Next, the child branches in the first branching order are attached to the main stem \(b_{\alpha}\), whose fitted centerline B-spline curve is given as \(\mathbf{C}^{\alpha}\). Let \(\nu_\delta\) be a junction node along the main stem path at which a path in \(\mathcal{G}_S\) branches off, and let \(\nu_\gamma\) be the distal node of that path; together, they define the child edge \(b_{\gamma}\). The node \(\nu_\delta\) corresponds to the parameter \(\bar{u}_\delta\) on the fitted B-spline curve \(\mathbf{C}^{\alpha}\). At two adjacent parameters on either side of \(\bar{u}_\delta\), \(\bar{u}_\delta^-\) and \(\bar{u}_\delta^+\), the knot multiplicity of the knot vector \(\mathcal{U}^{\alpha}\) is raised to \(p+1\) by repeated knot insertions. As a result, the B-spline curve \(\mathbf{C}^{\alpha}\) is split into three curves defined in the parametric intervals \([0, \bar{u}_\delta^-]\), \([\bar{u}_\delta^-, \bar{u}_\delta^+]\), and \([\bar{u}_\delta^+, 1]\). Because knot insertion is a geometry preserving operation [160], [162], the three curves together exactly reproduce the geometry of \(\mathbf{C}^{\alpha}\). The B-spline curve defined in the parametric interval \([\bar{u}_\delta^-, \bar{u}_\delta^+]\) is deleted to provide sufficient room for subsequent junction blending. This step serves a role similar to the junction pruning procedure carried out in Section 3.2.1. At initialization, the edge set in the branch graph \(\mathcal{G}_B\) is given as \(\mathcal{E}_B=\{b_{\alpha}\}\), and it is updated as follows. A junction node \(\nu_\delta\) is added to the main stem \(b_{\alpha}\). The edge \(b_{\alpha}\) is now defined between the nodes \(\nu_\delta\) and \(\nu_{\alpha}\), and its centerline B-spline curve is defined on the parameter interval \([\bar{u}_\delta^+, 1]\). A new edge, \(b_{\delta}\), is defined between the nodes \(\nu_{\mathrm{root}}\) and \(\nu_{\delta}\), with its corresponding B-spline curve defined on the parameter interval \([0, \bar{u}_\delta^-]\). A child edge \(b_{\gamma}\) is then defined from \(\nu_\delta\) to \(\nu_\gamma\) in the updated branch graph, with \(\mathcal{E}_B=\{b_{\delta}, b_{\alpha}, b_{\gamma}\}\).

The remaining branches are added similarly to the branch graph by carrying out the knot insertion and child branch attachment steps. Each edge with a newly fitted B-spline curve subsequently serves as the parent edge for its descendants, and the sub-paths at the junction nodes along this path are attached as its child edges. The path traversal proceeds in a depth-first manner to generate the complete branch graph. Thus, when a child edge is added, all the child edges of that child are processed first before moving to the next child of the same parent. Each child branch attachment operation adds a junction node, \(\nu_\delta\), creates a new edge, \(b_{\delta}\), updates the edge, \(b_{\alpha}\), and finally adds a child edge, \(b_{\gamma}\). Thus, the branch graph \(\mathcal{G}_B\) grows recursively from a single main stem edge into the complete branch graph, terminating once every point in the skeleton point cloud is fitted to the centerline B-spline curve. The advantage of the hierarchical construction strategy is that the knot insertion operation subdivides a B-spline curve without changing its geometry. Thus, regardless of how many times the parent B-spline curve has been subdivided to attach child branches, it still maintains its geometry and intrinsic smoothness while being represented as multiple B-spline curve segments compared to the global construction strategy. The hierarchical construction strategy is described in Algorithm 4.

3.3 B-spline volumetric parameterization for each branch↩︎

Each branch edge \(b_\alpha \in \mathcal{E}_{B}\) of the branch graph \(\mathcal{G}_B\) is represented by a cylindrical trivariate tensor-product B-spline volume \(\mathbf{V}^{\alpha}:[0,1]^3\rightarrow\mathbb{R}^3\) as \[\mathbf{V}^\alpha(r, \theta, u) \;=\; \sum_{i=1}^{\ell}\sum_{j=1}^{m}\sum_{k=1}^{n^\alpha} N_{i,s}(r)\, N_{j,q}(\theta)\, N_{k,p}(u)\, \mathbf{P}^\alpha_{i,j,k}, \quad (r, \theta, u)\in[0,1]^3, \label{eq:volume}\tag{5}\] where \(r\), \(\theta\), and \(u\) are the parametric coordinates, and \(s\), \(q\), and \(p\) are the degrees of the B-spline basis functions in the radial, circumferential, and axial directions, respectively. The B-spline control lattice \(\mathbf{P}^\alpha\) is of the size \((\ell \times m \times n^\alpha)\). \(\ell\) and \(m\) are kept constant for all branches to facilitate better junction blending between the incident volumes. Each basis function \(N_{i, s}(r)\), \(N_{j, q}(\theta)\), and \(N_{k, p}(u)\) is defined on the clamped knot vector \(\mathcal{R}\), the unclamped knot vector \(\Theta\), and the clamped knot vector \(\mathcal{U}^\alpha\), respectively. The circumferential knot vector, \(\Theta\), is unclamped in order to define periodic B-splines along the circumferential direction. To construct periodic B-splines, the last \(q\) columns are overlapped with the first \(q\) columns of the control lattice. Thus, \[\mathbf{P}^\alpha_{i,\, m - q + w,\, k} \;=\; \mathbf{P}^\alpha_{i,\, w,\, k} \;, \quad w = 1, 2, \dots, q, \label{eq:periodicity}\tag{6}\] ensures \(C^{q-1}\) surface continuity in the circumferential direction [136]. The tensor-product B-spline volume for each edge \(\mathbf{V}^\alpha\) is thus constructed by sweeping the \((r, \theta)\) parameterized cross-sections along the centerline curve \(\mathbf{C}^\alpha\) defined in Section 3.2 (Fig. 5 (a)).

A local orthonormal frame \(\{\mathbf{T}_k, \mathbf{N}_k, \mathbf{B}_k\}_{k=1}^{n^\alpha}\) is constructed along the centerline curve \(\mathbf{C}^\alpha\) at each Greville abscissa \(\{\xi_k\}_{k=1}^{n^\alpha}\) (Fig. 5 (b)), defined on \(\mathcal{U}^{\alpha}\) [160]. The unit tangent vector \(\mathbf{T}_k\) is tangent to the centerline curve at \(\mathbf{C}^{\alpha}(\xi_k)\), whereas the normal and binormal vectors, \(\mathbf{N}_k\) and \(\mathbf{B}_k\), span a cross-sectional plane that is perpendicular to \(\mathbf{T}_k\). The tangent vector is evaluated as \(\mathbf{T}_k = \dot{\mathbf{C}}^\alpha(\xi_k) / \|\dot{\mathbf{C}}^\alpha(\xi_k)\|\), where \(\dot{\mathbf{C}}^\alpha(\xi_k) =d\mathbf{C}^\alpha(\xi_k) / {du}\). The unit binormal vector is evaluated as the cross-product of the tangent and normal vectors, \(\mathbf{B}_k = \mathbf{T}_k \times \mathbf{N}_k\). The first unit normal vector along the centerline curve \(\mathbf{N}_1\) is obtained by projecting a fixed global up-vector, \((0, 0, 1)\), onto the plane perpendicular to \(\mathbf{T}_1\). If \(\mathbf{T}_1\) is nearly parallel to the default up-vector, an alternative up-vector is used. The remaining normals along the centerline curve are propagated using the double-reflection parallel transport method [163], producing a rotation-minimizing frame by minimizing twist about the tangential direction.

The control points of \(\mathbf{V}^\alpha\) are evaluated by sweeping circular cross-sections along the centerline curve \(\mathbf{C}^\alpha\) at the Greville abscissae \(\{\xi_k\}_{k=1}^{n^\alpha}\). At each \(\xi_k\), we define an angular location, \(\phi_j\), evaluated via uniform sampling over the circumference of the cross section. Similarly, at each \(\xi_k\), a radial location is defined as \(\rho_{i}\), evaluated through uniform sampling between the inner and outer radii of the cross section. This is given as \[\begin{align} \phi_j &= \frac{2\pi (j - 1)}{m - q}, \quad j = 1, \dots, m - q, \\ \rho_{i} &= h_k \left[ \eta + (1 - \eta) \tfrac{i-1}{\ell-1} \right], \quad i = 1, \dots, \ell, \label{eq:angular95radial95sampling} \end{align}\tag{7}\] where \(h_k\) denotes the outer branch radius at \(\xi_k\), obtained by interpolating the outer radius \(R_l\) at each skeleton node \(\nu_l\) along the centerline curve \(\mathbf{C}^\alpha\). The ratio of the inner branch radius to the outer branch radius, \(\eta \in (0, 1)\), is fixed for all branches. Each control point in the control lattice \(\mathbf{P}^{\alpha}\) is thus given as: \[\mathbf{P}^{\alpha}_{i,j,k} = \mathbf{C}^{\alpha}(\xi_k) + \rho_{i} \left( \cos(\phi_j)\,\mathbf{N}_k + \sin(\phi_j)\,\mathbf{B}_k \right), \qquad \begin{align} &i=1,\ldots,\ell,\\ &j=1,\ldots,m-q,\\ &k=1,\ldots,n^\alpha. \end{align} \label{eq:cp95construction}\tag{8}\] The \(q\) columns in the control-point lattice with circumferential index \(j = m - q + 1, \dots, m\) are evaluated according to the periodicity condition in Eq. 6 . By sweeping the cross-sections along the centerline and assembling the 3D control point lattice, we can then evaluate the corresponding tensor-product B-spline volume that produces the hollow cylindrical branch geometry shown in Fig. 5 (c-d).

3.4 Junction blending for continuity enforcement↩︎

As described in Section 3.3, each plant branch is represented by an independent B-spline volume constructed around its own centerline curve. Consequently, at a junction node \(\nu_\delta\in\mathcal{V}_B^{\mathrm{junc}}\), the junction-facing boundary cross-sections of the incident branch volumes are generally not coincident. The incident branch volumes must therefore be blended at the junction so that the assembled volumetric mesh contains no gaps and maintains smooth transitions across the junctions. In this section, we describe how these continuity conditions are enforced on the branch control lattices before mesh evaluation. We first describe the \(G^1\)-continuity requirements for a pair of B-spline curves in Section 3.4.1 and subsequently describe the extension for a pair of B-spline volumes in Section 3.4.2. We then introduce a cross-section splitting technique for the three incident branch control lattices meeting at bifurcation junctions in Section 3.4.3, and finally extend this approach to trifurcation junctions in Section 3.4.4. The continuity conditions enforced on the control lattices at each junction are reflected in the evaluated hex mesh. While this approach can be generalized to higher-valence junctions, extensions beyond trifurcations are outside the scope of the present work.

3.4.1 Junction continuity between two B-spline curves↩︎

Consider two clamped B-spline curves, \(\mathbf{C}^{\alpha}\) and \(\mathbf{C}^{\beta}\), shown in Fig. 6 (a), with control point sequences \(\mathbf{P}^\alpha_1, \dots, \mathbf{P}^\alpha_{n^{\alpha}}\) and \(\mathbf{P}^\beta_1, \dots, \mathbf{P}^\beta_{n^{\beta}}\), respectively. Here, \(n^{\alpha}\) and \(n^{\beta}\) denote the numbers of control points defining \(\mathbf{C}^{\alpha}\) and \(\mathbf{C}^{\beta}\), respectively. \(G^0\)-continuity is enforced by setting the terminal control point of \(\mathbf{C}^{\alpha}\) equal to the first control point of \(\mathbf{C}^{\beta}\). Thus, setting \(\mathbf{P}^\alpha_{n^{\alpha}} = \mathbf{P}^\beta_1\) closes the gap between the curves and enforces \(G^0\)-continuity, as shown in Fig. 6 (b). To achieve \(G^1\)-continuity, the tangent directions of the two curves must also be aligned at the junction. This can be imposed by making the shared control point and its immediate neighboring control points collinear [72], [160]. We impose both the continuity conditions simultaneously by placing the shared control point at the midpoint of the two adjacent control points: \[\mathbf{P}^\alpha_{n^{\alpha}} = \mathbf{P}^\beta_1 = \tfrac{1}{2}(\mathbf{P}^\alpha_{n^{\alpha}-1} + \mathbf{P}^\beta_2). \label{eq:g195midpoint}\tag{9}\] This ensures that the control points \(\mathbf{P}^\alpha_{n^{\alpha}-1}\), \(\mathbf{P}^\alpha_{n^{\alpha}}=\mathbf{P}^\beta_1\), and \(\mathbf{P}^\beta_2\) are collinear (Fig. 6 (c)), thereby providing a common tangent direction at the junction.

3.4.2 Junction continuity between two B-spline volumes↩︎

The \(G^1\)-continuity enforcement for curves in Eq. 9 extends naturally to a pair of adjacent trivariate tensor-product B-spline volumes. Consider two B-spline volumes, \(\mathbf{V}^{\alpha}\) and \(\mathbf{V}^{\beta}\), defined by the 3D control lattices \(\mathbf{P}^\alpha\) and \(\mathbf{P}^\beta\), respectively. For each fixed radial index, \(i\), and circumferential index, \(j\), the control points along the axial direction form an axial control column denoted by \(\mathbf{P}^{\alpha}_{i,j,:}\) or \(\mathbf{P}^{\beta}_{i,j,:}\), as shown in Fig. 7. The two B-spline volumes are joined by applying Eq. 9 column-wise to the axial control points for every fixed pair \((i, j)\). Specifically, \(G^0\)-continuity is enforced by making the terminal control point of each axial column in the control lattice of \(\mathbf{V}^{\alpha}\) equal to the first control point of each corresponding axial column in the control lattice of \(\mathbf{V}^{\beta}\). Thus, \(\mathbf{P}^{\alpha}_{i,j,n^{\alpha}} = \mathbf{P}^{\beta}_{i,j,1}\) for all \((i,j)\). This creates a shared boundary control cross-section between the two B-spline volumes, as shown in Fig. 7 (b). To enforce \(G^1\)-continuity, each control point on the shared boundary control cross-section is placed at the midpoint of its two neighboring axial control points, \(\mathbf{P}^\alpha_{i,j,n^{\alpha}-1}\) and \(\mathbf{P}^\beta_{i,j,2}\): \[\mathbf{P}^\alpha_{i,j,n^{\alpha}} = \mathbf{P}^\beta_{i,j,1} = \tfrac{1}{2}\bigl(\mathbf{P}^\alpha_{i,j,n^{\alpha}-1} +\mathbf{P}^\beta_{i,j,2}\bigr), \qquad \forall\,(i,j). \label{eq:g195volume95midpoint}\tag{10}\] Thus, three unique neighboring control points in each axial column become collinear, thereby providing a common axial tangent direction across the shared boundary between the two B-spline volumes.

3.4.3 Junction continuity at bifurcations↩︎

The pairwise continuity enforcement between two B-spline volumes described in Section 3.4.2 cannot be applied directly to bifurcations because three B-spline volumes meet at the junction. At a bifurcation junction node \(\nu_\delta \in \mathcal{V}_B^{\mathrm{junc}}\) in \(\mathcal{G}_B\), one parent edge \(b_\delta\) meets two child edges: \(b_{\gamma_1}\) and \(b_{\gamma_2}\). Let \(a \in \{\delta, \gamma_1, \gamma_2\}\) index the three branch edges incident at \(\nu_\delta\). For each incident edge \(b_a\), let \(\mathbf{C}^a\) denote its centerline curve, \(\mathbf{V}^a\) denote its tensor-product B-spline volume, and \(\mathbf{P}^a\) denote its corresponding control lattice. A control point in this lattice is denoted by \(\mathbf{P}^a_{i,j,k}\), where \(i\), \(j\), and \(k\) are the indices along the radial, circumferential, and axial directions, respectively. To enforce continuity at the bifurcation, the junction-facing control cross-section of each incident B-spline volume is first partitioned into two circumferential halves, as shown in Fig. 8. A splitting line is defined on each junction-facing control cross-section to partition the axial control columns into two circumferential halves. Corresponding axial control columns from paired circumferential halves are then blended using Eq. 10 . Because all incident control lattices have the same radial and circumferential resolution, corresponding circumferential halves can be paired across the three control lattices.

We define a junction plane that passes through the endpoints of the three centerline curves near the junction, as shown in Fig. 8 (a). Let \(\mathbf{p}^{\delta} = \mathbf{C}^{\delta}(1)\), \(\mathbf{p}^{\gamma_{1}} = \mathbf{C}^{\gamma_{1}}(0)\), and \(\mathbf{p}^{\gamma_{2}} = \mathbf{C}^{\gamma_{2}}(0)\) denote the junction-side endpoints of the parent and child curves, respectively. A unit normal vector to the junction plane, \(\hat{\mathbf{n}}\), is given as \[\hat{\mathbf{n}} \;=\; \frac{(\mathbf{p}^{\gamma_1} - \mathbf{p}^{\delta}) \times (\mathbf{p}^{\gamma_2} - \mathbf{p}^{\delta})}{\|(\mathbf{p}^{\gamma_1} - \mathbf{p}^{\delta}) \times (\mathbf{p}^{\gamma_2} - \mathbf{p}^{\delta})\|}. \label{eq:junction95vector}\tag{11}\] This unit normal vector provides a common reference direction for selecting the splitting line on each junction-facing control-point cross-section. For each incident branch, a control cross-section at a given axial index is assigned a depth index \(d\) based on its distance from the junction. Here, \(d=0\) denotes the junction-facing cross-section on each control lattice. For each incident edge \(a\), the axial index of each control point on the cross-section at depth-\(d\) is \[k^a_d = \begin{cases} n^a - d, & \text{if } a=\delta \\ 1 + d, & \text{if } a \in \{\gamma_1, \gamma_2\}, \end{cases} \label{eq:axial-depth-index}\tag{12}\] where \(n^a\) is the number of control points in the axial direction for the control lattice \(\mathbf{P}^{a}\). On the junction-facing control cross-section, with axial index \(k^{a}_{0}\), we define a radial line extending from \(\mathbf{P}^{a}_{1,j,k^{a}_{0}}\) to \(\mathbf{P}^{a}_{\ell,j,k^{a}_{0}}\). The splitting line is chosen as the radial diameter whose direction is most closely aligned with the normal vector \(\hat{\mathbf{n}}\). The corresponding circumferential index \(j^a_\mathrm{split}\) on the splitting line is thus given as: \[j^a_\mathrm{split} = \arg\max_{j\in\{1, \ldots, m-q\}} \left|\, \left(\frac{\mathbf{P}^{a}_{\ell,j, k^a_0} - \mathbf{P}^{a}_{1,j, k^a_0}}{\left\| \mathbf{P}^{a}_{\ell,j, k^a_0} - \mathbf{P}^{a}_{1,j, k^a_0} \right\|}\right)\cdot\hat{\mathbf{n}}\,\right|. \label{eq:splitting-index}\tag{13}\] The selected radial diameter partitions the control-lattice cross-section into two halves. Fig. 8 (c) shows the cross-sections partitioned by the splitting line on each control lattice cross-section, where matching colors indicate the corresponding pairs of axial columns used for pairwise continuity enforcement: \(b_{\delta}-b_{\gamma_1}\) (blue), \(b_{\delta}-b_{\gamma_2}\) (red), and \(b_{\gamma_1}-b_{\gamma_2}\) (green). For each paired set of axial columns, the continuity rule in Eq. 10 is applied column-wise. The control points on the depth \(d=0\) cross-section are placed at the midpoint of their corresponding depth \(d=1\) neighboring control points, thereby enforcing \(G^{1}\)-continuity across that interface.

The splitting line itself is treated as a common three-way control-lattice interface among the three incident branches at the junction. For each radial index \(i\) on the splitting line, we compute the centroid of the three adjacent depth \(d=1\) control points as \[\bar{\mathbf{P}}_i \;=\; \frac{1}{3} \sum_{a \in \{\delta, \gamma_1, \gamma_2\}} \mathbf{P}^{a}_{i, j^{a}_{\text{split}}, k^a_1}. \label{eq:junction-average}\tag{14}\] The control points on the splitting line of all three incident branches are then assigned to this centroid as \[\mathbf{P}^{a}_{i, j^{a}_{\text{split}}, k^a_0} = \bar{\mathbf{P}}_i, \qquad \forall\, a \in \{\delta, \gamma_1, \gamma_2\}. \label{eq:junction-assign}\tag{15}\] This operation assigns a common column of control points to the shared splitting line among the parent and child branch control lattices while using the adjacent depth \(d=1\) control columns to define the local tangent alignment at the junction.

3.4.4 Junction continuity at trifurcations↩︎

The enforcement of junction continuity at bifurcations, as described in Section 3.4.3, extends naturally to trifurcation junctions. At a trifurcation junction node \(\nu_\delta \in \mathcal{V}_B^{\mathrm{junc}}\) in \(\mathcal{G}_B\), one parent edge \(b_\delta\) meets three child edges \(b_{\gamma_1}\), \(b_{\gamma_2}\), and \(b_{\gamma_3}\). Let \(a \in \{\delta, \gamma_1, \gamma_2, \gamma_3\}\) index the four branch edges incident at \(\nu_\delta\). For each incident edge \(b_a\), let \(\mathbf{C}^a\) denote its B-spline centerline curve, \(\mathbf{V}^a\) denote its B-spline volume, and \(\mathbf{P}^a\) denote its corresponding control lattice. As in the case of bifurcation junctions, the junction-facing cross-section of each incident control lattice is split into two halves by selecting a radial diameter, also known as the splitting line. The circumferential index of the splitting line, \(j_{\mathrm{split}}^a\), is determined using Eq. 13 . The resulting partitioned control-lattice cross-sections are paired across the four incident edges as follows: \(b_{\delta}-b_{\gamma_2}\), \(b_{\gamma_2}-b_{\gamma_1}\), \(b_{\gamma_1}-b_{\gamma_3}\), and \(b_{\gamma_3}-b_{\delta}\), as shown in Fig. 9. For each paired set of axial columns in the control lattices, continuity is enforced using Eq. 10 .

The continuity enforcement procedure at the trifurcation junction differs from that at the bifurcation junction in two ways: the construction of the unit normal vector \(\hat{\mathbf{n}}\) and the treatment of the shared splitting-line interface. At a bifurcation junction, the three junction-side endpoints of the centerline curves define a unique junction plane. However, at a trifurcation junction, the endpoints \(\mathbf{p}^a\) of the four centerline curves \(\mathbf{C}^a\) are not always coplanar. Therefore, \(\hat{\mathbf{n}}\) is computed as the unit normal vector to the best-fit plane passing through these four endpoints. Let \[\mathbf{p}^{\delta} = \mathbf{C}^{\delta}(1), \qquad \mathbf{p}^{\gamma_r} = \mathbf{C}^{\gamma_r}(0), \quad r=1,2,3,\] denote the junction-side endpoints of the parent and child centerline curves. The centroid of these four endpoints is \(\bar{\mathbf{p}} = \frac{1}{4} \sum_{a \in \{\delta, \gamma_1, \gamma_2, \gamma_3\}} \mathbf{p}^{a}\). We define a \(4 \times 3\) matrix \(\mathbf{E}\) whose rows are \((\mathbf{p}^{a} - \bar{\mathbf{p}})^{\top}\) for \(a \in \{\delta, \gamma_1, \gamma_2, \gamma_3\}\). Applying singular value decomposition (SVD) to \(\mathbf{E}\), we obtain \(\hat{\mathbf{n}}\) as the right singular vector corresponding to the smallest singular value. The circumferential index on the splitting line \(j^a_\mathrm{split}\) is then determined for each control-lattice cross-section using Eq. 13 .

At a trifurcation junction, the splitting line is treated as a common four-way control-lattice interface among the incident B-spline volumes. For each radial index \(i\) on the splitting line, we compute the centroid of the four adjacent depth \(d=1\) control points as \[\bar{\mathbf{P}}_i = \frac{1}{4}\sum_{a\in\{\delta,\gamma_1,\gamma_2,\gamma_3\}}\mathbf{P}^{a}_{i,\,j^{a}_\mathrm{split},\,k^{a}_1}. \label{eq:centroid}\tag{16}\] The control points on the shared splitting line of all four incident branches are then assigned to this centroid: \[\mathbf{P}^a_{i, j^a_\mathrm{split}, k^a_0} = \bar{\mathbf{P}}_i, \quad \forall a \in \{\delta, \gamma_1, \gamma_2, \gamma_3\}. \label{eq:trifurcation-assign}\tag{17}\] However, the four adjacent depth \(d=1\) splitting-line control points, \(\mathbf{P}^{a}_{i,j^a_\mathrm{split},k^a_1}\), where \(a \in \{\delta, \gamma_1, \gamma_2, \gamma_3\}\), need not be coplanar. For each radial index \(i\), we construct a \(4 \times 3\) centered matrix whose rows are the offsets of these four points from their centroid: \[\mathbf{M}_i = \begin{pmatrix} \mathbf{P}^\delta_{i, j^\delta_\mathrm{split}, k^\delta_1} - \bar{\mathbf{P}}_i \\ \mathbf{P}^{\gamma_1}_{i,\,j^{\gamma_1}_\mathrm{split},\,k^{\gamma_1}_1} - \bar{\mathbf{P}}_i \\ \mathbf{P}^{\gamma_2}_{i,\,j^{\gamma_2}_\mathrm{split},\,k^{\gamma_2}_1} - \bar{\mathbf{P}}_i \\ \mathbf{P}^{\gamma_3}_{i,\,j^{\gamma_3}_\mathrm{split},\,k^{\gamma_3}_1} - \bar{\mathbf{P}}_i \end{pmatrix} \label{eq:adjacent95points95svd}\tag{18}\] Applying SVD to \(\mathbf{M}_i\), we obtain the local unit normal vector \(\hat{\mathbf{m}}_i\) as the right singular vector corresponding to the smallest singular value. The four adjacent depth \(d=1\) control points are then projected onto this best-fit plane. We first evaluate \(\tilde{\mathbf{P}}^{a}_{i}\) as \[\tilde{\mathbf{P}}^{a}_{i} = \mathbf{P}^{a}_{i,\,j^a_\mathrm{split},\,k^a_1} - \Bigl[\hat{\mathbf{m}}_i \cdot \Bigl(\mathbf{P}^{a}_{i,\,j^a_\mathrm{split},\,k^a_1} - \bar{\mathbf{P}}_i\Bigr)\Bigr]\hat{\mathbf{m}}_i, \qquad a \in \{\delta, \gamma_1, \gamma_2, \gamma_3\}. \label{eq:assignment}\tag{19}\] We then replace the original depth \(d=1\) control point \(\mathbf{P}^{a}_{i,\,j^a_\mathrm{split},\,k^a_1}\) with its projected point \(\mathbf{P}^a_{i, j^a_\mathrm{split}, k^a_1} = \tilde{\mathbf{P}}^{a}_{i}\). This coplanarity condition serves as the multi-way analog of the collinearity condition in the pairwise case (Eq. 10 ), and enables common tangent-plane alignment across the four incident B-spline volumes at the junction. Eqs. 16 –19 assign a common column of control points to the shared splitting line of the parent and child control lattices while enforcing coplanarity with the adjacent depth \(d=1\) control points that determine the local tangent alignment at the junction interface.

3.4.5 Hexahedral mesh generation↩︎

To generate the hex mesh vertices, each B-spline volume \(\mathbf{V}^\alpha\) is evaluated on a structured 3D parametric grid following Eq. 5 . The parametric grid is of the size \(S_r \times S_\theta \times S^\alpha_u\), where \(S_r\), \(S_\theta\), and \(S^\alpha_u\) denote the number of parametric coordinates along the radial, circumferential, and axial directions, respectively. The B-spline volume \(\mathbf{V}^\alpha\) of each branch edge is then evaluated at these parametric coordinates to generate a structured grid of mesh vertices. By construction, the control points in the control lattice of each B-spline volume are obtained by uniform sampling along the radial and circumferential directions (Eqs. 7 8 ). Therefore, the uniform sampling of radial and circumferential parametric coordinates results in uniformly placed mesh vertices in each branch cross-section. \(S_r\) and \(S_\theta\) thus remain constant for all the branch edges, which is important for enforcing junction continuity. However, along the axial direction, the control points are not uniformly distributed across all the branch edges. Therefore, the parametric coordinates in the axial direction are chosen so that mesh vertices have approximately equal spacing along the branch centerline curve. The number of parametric coordinates along the axial direction \(S^\alpha_u\) varies according to the branch length \(L^\alpha\). For a fixed spacing between the mesh vertices along the axial direction, \(\Delta_u\), we calculate \(S^\alpha_u\) as \[S^\alpha_u = \max\left(S_{\min}, \left\lceil \frac{L^\alpha}{\Delta_u} \right\rceil\right),\] where \(S_{\min}\) is a prescribed lower bound for the number of axial parametric coordinates, ensuring that short branches still contain a minimum number of mesh vertices along the axial direction. The control lattices, after enforcing junction blending for bifurcations and trifurcations, are shown in Fig. 9 (c), and the corresponding evaluated B-spline volumes are shown in Fig. 9 (d). In the circumferential direction, periodic B-splines defined on an unclamped knot vector are non-interpolatory at the control points along the shared splitting-line interface between the incident branches. Therefore, a small residual gap may remain between the evaluated B-spline volumes along the splitting-line interface. This gap is closed during hex mesh construction by averaging corresponding adjacent mesh vertices, producing a conforming hex mesh across each junction, as shown in Fig. 9 (e).

3.5 Mesh quality study at junctions↩︎

We evaluate the performance of the junction continuity enforcement by varying the angles between the incident branches for bifurcation and trifurcation junctions. For each generated hex mesh, we report the element quality using the scaled Jacobian metric. For an element \(e\), let \(\mathbf{t}_{A,1}\), \(\mathbf{t}_{A,2}\), and \(\mathbf{t}_{A,3}\) denote the three edge vectors meeting at a corner \(A\), ordered consistently with the right-handed reference hexahedron. Let \(L_{A,1}, L_{A,2}, L_{A,3}\) denote the corresponding edge lengths. The Jacobian determinant at the corner \(A\) is the signed volume of the parallelepiped spanned by these edge vectors, \(J_A = \det[\mathbf{t}_{A,1}, \mathbf{t}_{A,2}, \mathbf{t}_{A,3}]\). The scaled Jacobian at a corner \(A\) of the element is evaluated by normalizing \(J_A\) by the product of the three incident edge vectors. The scaled Jacobian of the element \(e\) is thus given as \[q_{e} = \min_{A=1,\ldots,8} \frac{J_A}{L_{A,1}L_{A,2}L_{A,3}}. \label{eq:min95scaled95jacobian}\tag{20}\] This metric satisfies \(|q_{e}| \le 1\), with \(q_{e} = 1\) corresponding to an orthogonal, positively oriented element. \(q_{e} \le 0\) indicates that at least one corner of the element is degenerate or inverted.

Fig. 10 shows different bifurcation junction configurations generated by varying the angle between the three incident branches. For each configuration, the angle between each branch is set to be between \(30°\) and \(270°\). For all the bifurcation junction configurations, we report positive scaled Jacobian values, as seen in Table 1. In addition, the \(10^\text{th}\) percentile of the element quality (\(P_{10}\)) satisfies \(P_{10} \ge 0.925\), indicating that at least \(90\%\) of the mesh elements in each configuration have \(q_{e} \ge 0.925\). Similarly, we perform the angle parametric study for trifurcation junctions, as shown in Fig. 11, varying the angles between the four incident branches at the junction. For all the configurations, we report good mesh quality, where all the elements are non-inverted (\(q_{e}>0\)). We maintain \(P_{10} \ge 0.905\) for trifurcation junctions, indicating that at least \(90\%\) of the mesh elements in each configuration have \(q_{e} \ge 0.905\). In all the junction configurations, the lowest \(q_{e}\) value is observed near the junction. Particularly, in one of the trifurcation junction configurations, the minimum \(q_{e}\) value is observed to be \(0.165\), as shown in Table 1. This is primarily due to increased localized curvature and extreme distortion near the junction introduced when the angles between the branches become very small.

4pt

4 Results and Discussion↩︎

We evaluate the proposed procedural modeling framework with respect to volumetric mesh generation from point cloud data, dynamic mesh updating, and suitability for FEA. We test the framework on three different plant species with increasing geometric and topological complexity. We first demonstrate volumetric mesh generation from skeleton point clouds of mung bean (Section 4.1.1) and tomato (Section 4.1.2) plants. In Section 4.2, we demonstrate the performance of the framework for unstructured point clouds of walnut tree branches. Here, in addition to generating volumetric meshes from plant skeletons, we also estimate the branch radii from unstructured point clouds, enabling mesh generation for complex branching architectures with spatially varying radii. We then demonstrate the main advantage of the procedural framework over static mesh generation approaches: automated modeling and dynamic mesh updating during plant growth without reconstructing the entire mesh from scratch at each growth stage (Section 4.3). Finally, in Section 4.4, we solve a steady-state diffusion problem on the generated meshes to demonstrate their direct suitability for FEA.

4.1 Volumetric mesh generation from plant skeleton datasets↩︎

We demonstrate the procedural modeling framework by generating volumetric meshes from skeleton point clouds of mung bean and tomato plants. Since these input point clouds approximate only the centerline geometry of the plants and do not provide branch-radius information, we model each branch using a circular cross-section with a prescribed constant radius along its length. Specifically, we assign a constant outer branch radius \(R_{b}\) to every node \(\nu_l\in \mathcal{V}_S\) in the skeleton graph \(\mathcal{G}_{S}\), such that \(R_l=R_{b}\), and use a constant inner-to-outer branch radius ratio \(\eta\).

4.1.1 Mung bean plants (Vigna radiata)↩︎

Point clouds of several mung bean varieties were obtained via terrestrial laser scanning [50], and skeletons defining the plant medial axes were extracted from these point clouds using the method described in [164]. In all mung bean specimens, the resulting skeleton graphs contain only bifurcation junctions. Because the skeleton point clouds are relatively sparse, additional nodes are inserted along the skeleton graph edges by interpolation, so that the spacing between adjacent nodes does not exceed \(1.0\) mm. A uniform outer branch radius \(R_{b} = 1.5\) mm is assigned to all skeleton nodes, and \(\eta = 0.5\) provides an inner branch radius of \(0.75\) mm. In the volumetric parameterization for each branch, cubic B-spline basis functions are used in all three parametric directions. In the circumferential direction, we use periodic B-splines with \(14\) unique control points, whereas in the radial direction, we use \(4\) control points. The hex mesh vertices for each branch are evaluated on a parametric grid of \(S_\theta = 28\) and \(S_r = 4\) coordinates, with an axial parametric spacing of \(\Delta_u = 1\,\text{mm}\). We use the hierarchical construction strategy to preserve the geometric smoothness of the parent centerline B-spline curves during branch insertion, as outlined in Algorithm 4. For each branch edge \(b_{\alpha}\), \(n^{\alpha}\) is set equal to the number of coordinates in \(\mathbf{Q}^{\alpha}\), as described in Section 3.2.2. We use the same parameters in the procedural framework for all the mung bean morphologies in the dataset, ranging from upright stem architectures to semi-erect, laterally sprawling branching architectures.

Fig. 12 (a) shows the volumetric meshes generated from the mung bean plant skeletons, along with the corresponding close-up views of selected bifurcation junctions colored by the element-wise scaled Jacobian \(q_{e}\) (Fig. 12 (b)). The minimum, mean, and \(10^{th}\) percentile (\(P_{10}\)) values of \(q_{e}\) are reported in Table ¿tbl:tab:mungbean? for each specimen. For all meshes, the mean \(q_{e}\) value is greater than \(0.97\), and the \(P_{10}\) value is greater than \(0.95\), indicating that at least \(90\%\) of the mesh elements have \(q_{e} > 0.95\). The minimum \(q_{e}\) values are approximately between \(0.1\) and \(0.27\), and the lower-quality elements are localized to a small number of junctions. These junctions typically have child branches meeting the parent branch at very acute angles. The close-up views of these junctions show visually smooth parent-child branch interfaces without any kinks on the surface, consistent with the continuity conditions enforced during junction blending.

4.1.2 Tomato plants (Solanum lycopersicum)↩︎

We next apply the procedural framework to four skeleton point clouds of tomato plants from the TomatoWUR dataset [165]. The selected specimens have dataset serial numbers 101, 336, 345, and 98, corresponding respectively to Tomato 1–4 in Fig. 13 and Table [tbl:tab:tomato]. Compared with the mung bean specimens, the tomato specimens contain more complex branching structures. In addition to plant stems, the skeletons also include large leaf-vein structures and smaller vein branches that attach near the distal ends of larger leaf branches. The knot insertion operation used in the hierarchical construction strategy (Section 3.2.2) requires nondegenerate parametric intervals on both sides of the child-branch attachment point on the parent B-spline curve. Specifically, if a child branch attaches to the parent B-spline centerline curve \(\mathbf{C}^{\alpha}(u)\) of branch \(b_{\alpha}\) at parameter \(\bar{u}_{\delta}\), the construction requires two adjacent parameters, \(\bar{u}_{\delta}^{-}\) and \(\bar{u}_{\delta}^{+}\), such that \(\bar{u}_{\delta}^{-}<\bar{u}_{\delta}<\bar{u}_{\delta}^{+}\). When the attachment point lies too close to the distal end of the parent B-spline curve, \(u=1\), there is insufficient parametric space to define a distinct neighboring parameter \(\bar{u}_{\delta}^{+}\) satisfying \(\bar{u}_{\delta}<\bar{u}_{\delta}^{+}<1\). Therefore, the hierarchical construction strategy is not suitable for this dataset because several smaller vein branches attach near the distal ends of their parent branches. We thus use the global construction strategy described in Section 3.2.1 to generate the B-spline centerline curves for the complete branch graph.

We increase the node density in each skeleton graph to achieve a maximum spacing of \(1.0\) mm between adjacent nodes. A constant outer branch radius \(R_{b} = 2\) mm is assigned to all the skeleton nodes, and \(\eta\) is set to a constant value of \(0.5\), resulting in an inner branch radius of \(1\) mm. The volumetric parameterization for each branch uses cubic B-spline basis functions in all three parametric directions. The number of unique control points in the circumferential direction is set to \(14\), and the number of control points in the radial direction is set to \(4\). The number of axial control points is determined from the branch length using the global construction strategy (Algorithm 3). The hex mesh vertices are evaluated with \(S_\theta = 14\) and \(S_r = 3\) parametric coordinates in each cross-section, and with an axial parametric spacing of \(\Delta_u = 2\,\text{mm}\). Fig. 13 shows the volumetric meshes for the tomato specimens, together with close-up views of selected junctions colored by the element-wise scaled Jacobian \(q_{e}\). As shown in Table [tbl:tab:tomato], the mean \(q_{e}\) value for each mesh is approximately \(0.95\), and \(P_{10}(q_{e})\) is approximately \(0.91\), both of which are slightly lower than the corresponding values obtained for the mung bean meshes. The reduction in mesh quality relative to the mung bean dataset is consistent with the higher branching density and greater variation in branch lengths of the tomato specimens. Moreover, a larger number of junctions increases the number of mesh elements affected by junction blending, which can locally reduce mesh quality. However, since each junction only affects a small number of elements, the overall mesh quality decreases very slightly.

4.2 Walnut tree mesh generation with spatially varying branch radii↩︎

The meshes generated for the mung bean and tomato plants in Section 4.1 use a constant outer branch radius for all branches within each specimen, independent of branch order. Real plants, particularly trees, exhibit substantial variation in branch radius across branch orders and along individual branches. Thus, the proposed framework must also support the blending of incident branches with different radii at junctions. We demonstrate this capability using three walnut tree specimens from the tree dataset [166], which provides unstructured 3D point clouds of branching structures together with their corresponding ground-truth skeletons. At each node \(\nu_l\) on the skeleton graph, we estimate the outer radius of cross-section \(R_{l}\) from the neighboring points in the unstructured point cloud. To compute \(R_l\), we first identify the points lying on the outer branch surface. A \(k\)-nearest-neighbor graph of the point cloud is constructed, and points with one-sided neighborhoods are classified as outer-surface points. The radius \(R_l\) at each node is then computed as the average distance from \(\nu_{l}\) to its \(k\) nearest outer-surface points. Finally, the radii \(R_{l}\) are interpolated along the fitted centerline B-spline curve \(\mathbf{C}^{\alpha}(u)\) of each branch to obtain the outer branch radius \(h_{k}\) at each Greville abscissa \(\xi_{k}\) when constructing the control lattice in Eq. 8 . In this way, the spatially varying branch radius along each centerline B-spline curve is incorporated directly into the volumetric parameterization, without requiring control-point optimization. The sparse plant skeletons provided in the dataset are resampled to a maximum internodal spacing of \(1\,\text{mm}\), yielding a denser skeleton point cloud for more accurate B-spline centerline curve fitting. The centerline curves for the complete branch graph are constructed using the hierarchical construction strategy described in Section 3.2.2, preserving the geometric smoothness of the parent centerline B-spline curve each time a child branch is attached to it. The volumetric parameterization uses cubic basis functions in all three parametric directions, with \(14\) unique control points in the circumferential direction and \(4\) control points in the radial direction. The hex mesh vertices are computed by evaluating each B-spline volume on a cross-sectional parametric grid of size \(S_\theta = 42\) and \(S_r = 4\), with the axial parametric spacing \(\Delta_u\) set per specimen based on its branch-length distribution.

The reconstructed volumetric meshes of the walnut tree specimens are shown in Fig. 14 (a), with close-up views of representative junction regions in Fig. 14 (b). The input unstructured point cloud is overlaid on each mesh to visualize reconstruction accuracy, with each point colored by its reconstruction error, defined as the shortest Euclidean distance from the point to the mesh. Table [tbl:tab:walnut95reconstruction95error] summarizes the reconstruction error statistics across all tree specimens. We report a mean error value of less than \(3\) mm and a maximum error value not exceeding \(38.6\) mm across all the specimens. As reported in Table ¿tbl:tab:walnut?, the meshes for walnut tree specimens achieve the highest overall mesh quality among the three datasets, with mean \(q_{e}\) and \(P_{10}(q_{e})\) both greater than \(0.99\). This improvement relative to the meshes obtained in Section 4.1 is due to the finer circumferential discretization (\(S_\theta = 42\) for walnut, versus \(28\) for mung bean and \(14\) for tomato specimens), together with a larger number of mesh elements per branch generated in the taller walnut trees. Because the mesh elements with lower quality remain confined to a localized region at each junction, increasing the number of elements per branch subsequently increases the proportion of high-quality elements away from the junctions and thus improves the overall mesh quality. The mesh elements corresponding to the minimum \(q_{e}\) values remain localized to only a small number of junction regions where branches meet at acute angles.

6pt

4.3 Dynamic plant growth modeling↩︎

Modeling plant growth and development is a dynamic process in which the meshes must robustly adapt to extending and bifurcating branches that continuously modify the underlying plant topology. The computational cost of static mesh generation pipelines can compound over multiple growth stages, thus becoming extremely prohibitive for the dynamic modeling of complex branching architectures. In our procedural framework, we first represent the plant centerline geometry as a graph and subsequently fit it to smooth B-spline centerline curves. A tensor-product B-spline volume is then constructed around the centerline curve of each branch and finally assembled to generate a hex mesh after junction-blending operations. This sequence of steps allows plant growth to be modeled using local graph updates. Newly introduced branch segments are first converted to B-spline branch volumes, and junction blending is carried out in a local manner while retaining the hex mesh generated in the previous growth stage. This provides a scalable approach for dynamic volumetric modeling of growing plants, enabling local mesh updates without constructing geometry from scratch.

We demonstrate dynamic mesh updating using the Mung bean 3 specimen (Section 4.1.1) by simulating growth through a breadth-first search (BFS) traversal of the skeleton graph, starting from the root node and moving systematically according to the branching order. At each growth stage, only the part of the skeleton that has grown so far, referred to as the active skeleton subgraph, is considered for volumetric mesh generation. The active skeleton subgraph contains all nodes whose BFS distance from the root node is less than or equal to the current growth level. Fig. 15 shows the generated volumetric hex meshes at six representative growth stages. The hex mesh generated at previous growth stages is shown in blue, whereas newly added branches at the current growth stage are shown in orange. As reported in Table ¿tbl:tab:growth95study?, the hex meshes generated at all the growth stages have positive scaled Jacobian values, indicating that there are no inverted mesh elements. The mean \(q_{e}\) value remains at least \(0.9729\) across all growth stages, indicating that the local mesh updates at each stage preserve global mesh quality. Table [tbl:tab:growth95study95times] compares the computation time required to rebuild the full plant geometry at each growth stage with the time required for dynamic updating by locally adding the new branches. We demonstrate that reconstructing the entire plant geometry from scratch requires \(0.096\) seconds at the final growth stage (Step 5). In contrast, the dynamic approach only requires between \(0.03\) and \(0.05\) seconds to generate the volumetric mesh at each growth stage after the initial step (Step 0). This confirms that the framework decouples local growth updates from global branching complexity, providing a scalable framework for dynamic plant modeling.

4.4 Diffusion simulation in volumetric meshes↩︎

To evaluate the analysis-suitability of the generated volumetric meshes, we solved a steady-state diffusion equation on each plant geometry. The purpose of these simulations is to demonstrate the direct solver-readiness of the volumetric meshes generated by the proposed procedural modeling framework. We thus verify that the generated hex meshes are conforming and have valid element orientations. This is important for downstream simulations of plant biomechanics and transport phenomena with volumetric meshes, where poor mesh quality could affect the reliability of FEA.

We define a steady-state nutrient concentration field \(c(\mathbf{x})\) over the 3D plant domain \(\Omega \in \mathbb{R}^3\). We denote the unique branch incident to the root node \(\nu_{\mathrm{root}}\), which is the proximal node in the branch graph \(\mathcal{G}_{B}\), as \(b_{\zeta} \in \mathcal{E}_B\). The B-spline volume \(\mathbf{V}^{\zeta}\) is defined for the branch \(b_{\zeta}\). A concentration \(c(\mathbf{x})=1\) is prescribed on the root (source) cross-section \(\Gamma_{\mathrm{root}}\), defined as the proximal (\(u=0\)) face of \(\mathbf{V}^{\zeta}\): \[\Gamma_{\mathrm{root}} = \bigl\{\, \mathbf{V}^{\zeta}(r,\theta,0) \;:\; r\in[0,1],\;\theta\in[0,1] \,\bigr\}.\] This cross-section is referred to as the source boundary cross-section. A constant outward flux is prescribed on the terminal leaf (sink) cross-sections, defined as the distal (\(u=1\)) faces of all leaf branches. Let \(\mathcal{I} = \bigl\{\, \alpha \;:\; b_\alpha \in \mathcal{E}_B,\; \nu_\alpha \in \mathcal{V}_B^{\mathrm{leaf}} \,\bigr\}\) denote the index set of leaf branches. We then define the terminal leaf cross-sections as \[\Gamma_{\mathrm{leaf}} = \bigcup_{\alpha \in \mathcal{I}} \bigl\{\, \mathbf{V}^{\alpha}(r,\theta,1) \;:\; r\in[0,1],\;\theta\in[0,1] \,\bigr\}.\] The lateral branch walls are the inner and outer cylindrical surfaces of every branch volume, which are the faces at \(r=0\) (inner) and \(r=1\) (outer) of each \(\mathbf{V}^{\alpha}\): \[\Gamma_{\mathrm{wall}} = \bigcup_{b_\alpha \in \mathcal{E}_B} \Bigl( \bigl\{\, \mathbf{V}^{\alpha}(0,\theta,u) \;:\; \theta\in[0,1],\;u\in[0,1] \,\bigr\} \;\cup\; \bigl\{\, \mathbf{V}^{\alpha}(1,\theta,u) \;:\; \theta\in[0,1],\;u\in[0,1] \,\bigr\} \Bigr).\] The governing problem is modeled as the Laplace equation, with a Dirichlet boundary condition at \(\Gamma_{\mathrm{root}}\) and an inhomogeneous Neumann (flux) boundary condition at \(\Gamma_{\mathrm{leaf}}\). A homogeneous Neumann (no flux) boundary condition is applied at \(\Gamma_{\mathrm{wall}}\). The diffusion problem is thus given as \[\begin{align} \nabla^2 c(\mathbf{x}) &= 0 && \text{in } \Omega, \\ c(\mathbf{x}) &= 1 && \text{on } \Gamma_{\mathrm{root}}, \\ \nabla c(\mathbf{x})\cdot\mathbf{n} &= -100 && \text{on } \Gamma_{\mathrm{leaf}}, \\ \nabla c(\mathbf{x})\cdot\mathbf{n} &= 0 && \text{on } \Gamma_{\mathrm{wall}}, \end{align}\] where \(\mathbf{n}\) is the outward unit normal defined on the domain boundary \(\partial \Omega= \Gamma_{\mathrm{root}} \,\cup\, \Gamma_{\mathrm{leaf}} \,\cup\, \Gamma_{\mathrm{wall}}\). Although the steady-state diffusion problem does not completely capture the complexity of the nutrient transport process, it provides a useful test for ensuring the smooth propagation of the concentration solution field through multiple branching levels and across the blended junction regions.

We perform diffusion simulations on the hex meshes generated from three mung bean, tomato, and walnut tree specimens (Figs. 16– 18). The generated meshes span different levels of branching complexity, including meshes with uniform and spatially varying radii. The solver converges for every mesh, indicating that the meshes are valid for physics-based simulations. Across all specimens, the concentration field \(c(\mathbf{x})\) decreases monotonically from the root cross section toward the terminal leaf cross-sections and remains continuous across all branch junction regions. These results demonstrate the capability of the proposed procedural volumetric modeling framework to not only generate valid meshes but also to support their direct use in FEA. The solver-readiness of the volumetric modeling framework highlights its advantage over surface reconstruction methods, which often require additional mesh generation and cleanup pipelines to be carried out before physics-based simulations can be performed.

5 Conclusion and Future Work↩︎

We present a novel procedural framework that converts plant skeleton point clouds into analysis-suitable hexahedral meshes for plant branching structures. The pipeline extracts a branch graph from the skeleton, fits a B-spline centerline curve to each branch, constructs a cylindrical tensor-product B-spline volume along each centerline curve, and enforces smooth junction continuity across incident branch volumes at each junction via control point blending operations to generate the hex meshes. Mesh generation at several angle configurations for bifurcation and trifurcation junctions shows that mesh quality deteriorates only at extremely acute angles while remaining free of element inversions and intersections in all cases. In highly complex branching architectures, localized element inversions can occur near the junctions. However, these can be resolved by increasing the number of control points or by increasing the pruning radius near junctions to allow more room for junction blending. We demonstrate the generalizability of the framework across different plant species by generating analysis-suitable plant geometries with consistently good mesh quality. We also demonstrate that plant growth can be modeled dynamically by extending the skeleton graph while adding new branches at each growth stage and subsequently carrying out blending at the newly added junctions. The resulting meshes are solver-ready, as evidenced by a steady-state diffusion simulation in which the concentration field remains smooth across every junction for all the generated plant geometries. Together, these results demonstrate that the framework bridges the gap between surface-based plant reconstruction and analysis-ready procedural volumetric mesh generation pipelines. Several directions remain for future work. Our framework currently focuses on bifurcation and trifurcation junctions; however, some plant architectures exhibit higher-valence junctions. Extending the junction blending operations to such junctions would enable the modeling of more complex branching architectures. Moreover, while the present work uses B-spline volumes as a scaffold for generating conforming hexahedral meshes, we plan to develop a complete B-spline-based volumetric parameterization of the entire plant geometry for isogeometric analysis.

Acknowledgements↩︎

A. Moola and A. Pawar would like to acknowledge funding in part from the National Science Foundation under award number DMS-2533961 and the Translational AI Center at Iowa State University.

References↩︎