New Coupled Canopy–Light Model (CCLM) to Improve Visual Polymorphism Simulation of Fir Morphology

: Environmental factors substantially inﬂuence the growth of trees. The current studies on tree growth simulation have mainly focused on the effect of environmental factors on diameter at breast height and tree height. However, the inﬂuence of environmental factors, especially light


Introduction
Forest tree canopy shapes are diverse, and many factors contribute to their diversity.In addition to the genetic factors of the forest tree itself, the growing environment of the forest tree is also crucial for canopy shape diversity [1,2].Many environmental factors affect the canopy morphology of forest trees, including fire, water, light, temperature, soil, elevation, spatial structure, species interaction, and so on [3,4].At the scale of the sample plot, differences in soil, temperature, and elevation within the same sample plot are not significant [5][6][7].However, competition for spatial resources in forest stands can lead to significant differences in spatial structure between forest trees and their neighboring trees, resulting in significant differences in light conditions for forest trees within forest stands [8][9][10].The growth of the canopy varies under different light conditions.
Therefore, it is important to consider the effect of light on the canopy shape, and there are some scholars who have already started such research [11,12].Shanin et al. proposed a new model, which operates with the 3D-representation of tree canopies and light transmission through the canopy using discrete spatial and temporal resolution, paying special attention to the simulation of asymmetry in canopy shape as a result of competition for light [13].Duchemin et al. considered the growth of the tree canopy as a continuous front propagation process and propose a model which only requires two parameters to describe the shape of the crown: the intensities of phototropic and gravitropic growing responses [14].Hu et al. found that with the increase of light, the canopy of seedlings developed from wide and loose to relatively tight and narrow [15].He et al. proposed that shading can increase the leaf area of ginkgo seedlings, prolong their growth period, and promote the growth of long branches and plant height [16].Xu et al. believed that hard light promoted the differentiation of branches by inhibiting the growth of the trunk, and that shading promoted the rebranching of branches by inhibiting the growth of first-order branches [17].The canopy shape of a forest tree is closely related to its ability to obtain light resources, which affects its light distribution.The light distribution of the canopy also affects the canopy shape of the forest tree [18], so a mutual influence exists between the two.When the light intensity is below or above a certain value, branch growth stagnation or even death may occur.Between these limits, the growth rate of branches first increases and then decreases with the increase in light intensity [19][20][21].
To research the influence of light on the canopy shape, it is necessary to calculate the canopy light distribution [22][23][24].The methods used for calculating plant canopy light distribution are constantly advancing with the development of technology, from the direct observation of plants to measuring plants with the aid of instruments and, more recently, to the digital simulation of plants [25][26][27].In the process, the accuracy and efficiency of model calculations are constantly increased [28][29][30].Initially, canopy analyzers were used to capture photos of plants in different orientations, and physical indicators such as light spots and light transmission were analyzed based on the acquired images [31].Alternatively, a light meter can be used to determine the light distribution of a plant canopy [32], a method that relies on hardware devices to measure the light distribution information of plants, which is a cumbersome and inefficient process.Lao [33] and Wang [34] used a 3D digitizer to measure the 3D digital morphology of plant canopies and analyzed the canopy light distribution.However, 3D digitizers cannot rapidly acquire data and do not accurately capture the actual light distribution.
With the continuous development of information technology, hardware and software are more often being combined to construct mathematical physical models to study the light distribution in plant canopies and to construct 3D models of plants [35,36].Simulations are mainly performed using various illumination models.Illumination models, also known as light-dark effect models, are mathematical models that replace real physical models of light irradiation on the surface of an object.The two main types of illumination models are local [37] and global [38].Local illumination models can be subdivided into the Lambert, Gourand, and Phong illumination models [39,40].Global illumination models can be subdivided into the ray-tracing, radiance, and photon mapping algorithms [41][42][43].
Among them, the most widely studied is the ray-tracing algorithm.Zhang studied the key technology of the apple-tree canopy illumination model and implemented a GPU-based method to calculate plant canopy light [44].Hua studied tree canopy light distribution based on the 3D spatial structure and ray-tracing algorithm, scanned and reconstructed the 3D structure of poplars, and implemented the canopy light distribution algorithm based on the GPU [45].
In terms of studies on the effect of light on the canopy morphology of forest trees, Shanin et al. and Duchemin et al. quantitatively analyzed the effect of light on canopy morphology at the single-tree scale.Hu et al. qualitatively described light intensity using the words "forest gap", "full light", and "shade", and then qualitatively studied the effect of light on canopy morphology at the forest-stand scale.However, quantitative study of the mathematical relationship between light intensity and canopy width at the foreststand scale are lacking.Regarding plant canopy light distribution, the above-mentioned researchers have tended to focus on the canopy light distribution of a single tree, ignoring the light distribution in the whole forest stand.No further studies have been conducted on the relationship between light distribution and canopy morphology.
As the main management tree species of Huangfengqiao Forest Farm, fir trees have a wide range of ecological, economic and medicinal values [46].In this study, fir was selected as the research object.It is worth noting that the method in this paper also can be applied to other non-fir species.The purpose of our study was to put forward a new coupled canopylight model (CCLM) based on light distribution data and canopy morphological data.The CCLM expresses the mathematical relationship between light intensity and canopy width and was used to improve the visual polymorphism simulation of fir morphology.The development of the CCLM involved four steps: (1) constructing a 3D scene in UE4 by forestry survey data, (2) obtaining the light distribution data of the fir tree canopy in a 3D scene with a ray-tracing algorithm, (3) building a coupled canopy-light model (CCLM) of fir and performing a visual simulation, and (4) comparing the CCLM with the method proposed by Ma et al. [47] and calculating the accuracy of the CCLM in visual polymorphism simulation of fir morphology.

Study Area
We selected the Huangfengqiao Forestry Farm in China as the study area, which is in a zonal distribution located at the east and west border of You County in Hunan Province, with a latitude and longitude range of 113 • 04 to 113 • 43 E and 26 • 43 to 27 • 06 N. The terrain is dominated by low and medium mountainous landscapes, with slopes ranging from 20 to 35 degrees and elevation from 115 and 1270 m.The climate is subtropical monsoonal, with an average annual temperature of 17.8 • C, an annual average precipitation of 1410.8 mm, a frost-free period of 292 days, and an average annual sunshine time of 1612 h.The forest management area is 152,438 acres, with 891,262 cubic meters of forest wood storage.The forest is rich in trees and species, with more than 430 species of woody plants, including more than 10 species of rare and endangered trees under national protection, such as Taxus chinensis var.mairei, ginkgo, flame tree, and tulip tree [48] (Figure 1).

Data Sampling
Considering factors such as site conditions, altitude, slope, aspect, difficulty of data collection, etc., we selected a typical sample plot with an area of 40 m by 40 m and a total of 230 fir trees in the Huangfengqiao Forestry Farm.We used a UAV orthophoto to obtain the spatial topographic elevation of the sample plot.An electronic total station was used to obtain the relative coordinates (X, Y) of the roots of the fir trees to represent the distribution of the fir trees in the sample plot (Figure 2).The fir trees in the sample plot were surveyed for three consecutive years from 2015 to 2017 with conventional field surveys; leveling poles and laser altimeters were used to measure tree crown morphology data, and we recorded the tree height, canopy width, height at maximum canopy width, and under-living branch height in four quadrants (as shown in Table 1).The canopy width was calculated in one-meter intervals for the corresponding height in the four quadrants.We collected the global long-series high-resolution photosynthetically active radiation (PAR) data from the National Tibetan Plateau Science Data Center [49].

Data Sampling
Considering factors such as site conditions, altitude, slope, aspect, difficulty of data collection, etc., we selected a typical sample plot with an area of 40 m by 40 m and a total of 230 fir trees in the Huangfengqiao Forestry Farm.We used a UAV orthophoto to obtain the spatial topographic elevation of the sample plot.An electronic total station was used to obtain the relative coordinates (X, Y) of the roots of the fir trees to represent the distribution of the fir trees in the sample plot (Figure 2).The fir trees in the sample plot were surveyed for three consecutive years from 2015 to 2017 with conventional field surveys; leveling poles and laser altimeters were used to measure tree crown morphology data, and we recorded the tree height, canopy width, height at maximum canopy width, and under-living branch height in four quadrants (as shown in Table 1).The canopy width was calculated in one-meter intervals for the corresponding height in the four quadrants.We collected the global long-series high-resolution photosynthetically active radiation (PAR) data from the National Tibetan Plateau Science Data Center [49].

Data Sampling
Considering factors such as site conditions, altitude, slope, aspect, difficulty of data collection, etc., we selected a typical sample plot with an area of 40 m by 40 m and a total of 230 fir trees in the Huangfengqiao Forestry Farm.We used a UAV orthophoto to obtain the spatial topographic elevation of the sample plot.An electronic total station was used to obtain the relative coordinates (X, Y) of the roots of the fir trees to represent the distribution of the fir trees in the sample plot (Figure 2).The fir trees in the sample plot were surveyed for three consecutive years from 2015 to 2017 with conventional field surveys; leveling poles and laser altimeters were used to measure tree crown morphology data, and we recorded the tree height, canopy width, height at maximum canopy width, and under-living branch height in four quadrants (as shown in Table 1).The canopy width was calculated in one-meter intervals for the corresponding height in the four quadrants.We collected the global long-series high-resolution photosynthetically active radiation (PAR) data from the National Tibetan Plateau Science Data Center [49].

Canopy Shape Description and Light Parameters
The commonly used canopy shape descriptors include under-living branch height, height at maximum canopy width, canopy width, and so on.When using these indicators to identify similarities in canopy shapes, the accuracy is insufficient.We therefore addressed this issue by measuring the width of the canopy at one-meter intervals to obtain a more detailed description of the morphological canopy characteristics.When we compared the two groups of canopy morphology models, the canopy of each tree in the first group of models was recorded at one-meter intervals in the four quadrants to obtain a set of canopy widths x n .The canopy widths of the corresponding trees, as well as their directions and heights, were recorded for the second group of models to acquire a set of canopy widths y n (Figure 3).We compared the two sets of data x n and y n using the Euclidean distance to identify the similarity between the two groups in terms of the canopy shapes.The smaller the Euclidean distance, the more similar the canopy shapes between the two groups were.

Canopy Shape Description and Light Parameters
The commonly used canopy shape descriptors include under-living branch height, height at maximum canopy width, canopy width, and so on.When using these indicators to identify similarities in canopy shapes, the accuracy is insufficient.We therefore addressed this issue by measuring the width of the canopy at one-meter intervals to obtain a more detailed description of the morphological canopy characteristics.When we compared the two groups of canopy morphology models, the canopy of each tree in the first group of models was recorded at one-meter intervals in the four quadrants to obtain a set of canopy widths  .The canopy widths of the corresponding trees, as well as their directions and heights, were recorded for the second group of models to acquire a set of canopy widths  (Figure 3).We compared the two sets of data  and  using the Euclidean distance to identify the similarity between the two groups in terms of the canopy shapes.The smaller the Euclidean distance, the more similar the canopy shapes between the two groups were.There are two common methods of measuring light intensity: the radiant illuminance, which represents the radiant flux per unit area in watts per square meter (W/m 2 ), indicating the amount of energy received by the surface of the object; and the light There are two common methods of measuring light intensity: the radiant illuminance, which represents the radiant flux per unit area in watts per square meter (W/m 2 ), indicating the amount of energy received by the surface of the object; and the light intensity, which indicates the light flux per unit area in lux, indicating the amount of illumination of the surface of the object.Light intensity is the human visual perception of light, indicating light and darkness, whereas the perception of light by plants can be measured by radiation illumination.Solar radiation is not always effective for plant photosynthesis; the component of solar radiation that is effective for plant photosynthesis is called photosynthetically active radiation (PAR), which has a wavelength range of 380 to 710 nm.The annual average PAR is used as the light parameter in the ray-tracing algorithm.Because the light condition of a canopy is affected by the angle of solar irradiation, the angle of solar irradiation and PAR vary at different times.In a 3D scene, the changes in solar orientation and PAR are simplified by dividing them into seven moments (6:00, 8:00, 10:00, 12:00, 14:00, 16:00, and 18:00), and the sum of PAR at these seven moments at the canopy measurement point is calculated to represent the PAR at that point.

Global Illumination Algorithm Based on Ray Tracing
The ray-tracing algorithm is a global illumination algorithm that can accurately simulate global illumination.Unlike local illumination models, which only consider a single reflection or refraction from a light source to the human eye or camera, global illumination models account for the effect of other nonlight surfaces on the result.The ray-tracing algorithm simulates light in reality and tracks the path of light propagation.In a scene, after the light and object collide, a secondary light is generated due to reflection and refraction, and the secondary light enters the scene for collision testing until the number of iterations reaches the expected value, or the light is finally shot into the camera.Considering the probability of the light eventually entering the camera, combined with the principle of reversibility of the optical path, light is emitted from the camera to the light source and then backtracked to calculate the loss of energy, which is finally colored on a light screen.This operation is repeated for each pixel to obtain a 2D representation of a 3D scene.In this study, we used real survey data to construct a virtual forest stand environment, and we then used a ray-tracing algorithm to calculate the light distribution of each part of the canopy to provide light distribution data for the construction of the CCLM (Figure 4).intensity, which indicates the light flux per unit area in lux, indicating the amount of illumination of the surface of the object.Light intensity is the human visual perception of light, indicating light and darkness, whereas the perception of light by plants can be measured by radiation illumination.Solar radiation is not always effective for plant photosynthesis; the component of solar radiation that is effective for plant photosynthesis is called photosynthetically active radiation (PAR), which has a wavelength range of 380 to 710 nm.The annual average PAR is used as the light parameter in the ray-tracing algorithm.
Because the light condition of a canopy is affected by the angle of solar irradiation, the angle of solar irradiation and PAR vary at different times.In a 3D scene, the changes in solar orientation and PAR are simplified by dividing them into seven moments (6:00, 8:00, 10:00, 12:00, 14:00, 16:00, and 18:00), and the sum of PAR at these seven moments at the canopy measurement point is calculated to represent the PAR at that point.

Global Illumination Algorithm Based on Ray Tracing
The ray-tracing algorithm is a global illumination algorithm that can accurately simulate global illumination.Unlike local illumination models, which only consider a single reflection or refraction from a light source to the human eye or camera, global illumination models account for the effect of other nonlight surfaces on the result.The ray-tracing algorithm simulates light in reality and tracks the path of light propagation.In a scene, after the light and object collide, a secondary light is generated due to reflection and refraction, and the secondary light enters the scene for collision testing until the number of iterations reaches the expected value, or the light is finally shot into the camera.Considering the probability of the light eventually entering the camera, combined with the principle of reversibility of the optical path, light is emitted from the camera to the light source and then backtracked to calculate the loss of energy, which is finally colored on a light screen.This operation is repeated for each pixel to obtain a 2D representation of a 3D scene.In this study, we used real survey data to construct a virtual forest stand environment, and we then used a ray-tracing algorithm to calculate the light distribution of each part of the canopy to provide light distribution data for the construction of the CCLM (Figure 4).

Construction of Coupled Canopy-Light Model (CCLM)
The Coupled canopy-light model is shown in Figure 5. From the field surveys, we obtained base canopy morphological data (Base_CMDD), expressed as a set (i, d, h, w), which indicates the width of the canopy for the corresponding tree (number i), direction (direction d), and height (height h) as w.We used points (h, w) as the type value points (data points describing the geometry of the curve)  of the cubic uniform B-spline curve [50].The control points  were solved by the type value points.The system of equations to solve for the control points is shown in Equation (1).

Construction of Coupled Canopy-Light Model (CCLM)
The Coupled canopy-light model is shown in Figure 5. From the field surveys, we obtained base canopy morphological data (Base_CMDD), expressed as a set (i, d, h, w), which indicates the width of the canopy for the corresponding tree (number i), direction (direction d), and height (height h) as w.We used points (h, w) as the type value points (data points describing the geometry of the curve) Q i of the cubic uniform B-spline curve [50].The control points P i were solved by the type value points.The system of equations to solve for the control points is shown in Equation (1).The control points P i were counted based on the system of equations in (1), and the B-sample curve of the canopy curve was calculated using P i with the expression: For the cubic uniform B-spline curve, the basis function is F i,3 (t): The canopy curves of each tree in the four quadrants could be obtained, and the 3D model of the tree was constructed from the canopy curves in each quadrant, which we combined with the topographic data and the distribution information of the fir trees to construct the 3D scene of the sample plot.The light distribution data (LDD), expressed as (i, d, h, l), was calculated by combining the ray-tracing algorithm with the annual average PAR as the light parameter in the 3D scene.
The base period tree attribute data (Base_TAD) data are expressed as (i, D, H), where i represents the tree number, D represents the diameter at breast height, and H represents the tree height.We obtained the current tree attribute data, also expressed as (i, D, H), from a simulation (Current_TAD_Simulation) using the diameter in the breast-diameter growth model (BDGM) and the height-diameter curve model (HDCM).The BDGM and HDCM were based on a previous study [51], where RS is the relative plant distance, RD is the relative dominance, D is the diameter at breast height, and SI is the status index, reflecting the combination of topography, soil, and climate.HDCM : dD/dt = 1.564RS 0.515 RD 0.027 0.133SI 0.886 D 0.230 − 0.016SI 0.733 D (4) A quantitative relationship exists: from which the current canopy morphological data without the effect of light (Cur-rent_CMDD_N) is obtained, expressed as (i, d, h, w), where H denotes H in Base_TAD; H denotes H in Current_TAD_Simulation; and h and w denote h and w in Base_CMDD, respectively; h and w denote h and w in Current_CMDD_N, respectively.The current canopy morphological data (Current_CMDD) were obtained through field surveys, and the CCLM was finally constructed using statistical methods, expressed as the equation of Current_CMDD (i, d, h, w), Current_CMDD_N (i, d, h, w), and LDD (i, d, h, l).

Contrast Validation Model
After we constructed the CCLM, we incorporated the LDD and Current_CMDD_N into the CCLM to obtain the effects of light on the current canopy morphology, which is the current canopy morphological data obtained from the simulation using the CCLM (Current_CMDD_CCLM).
We then combined the data on under-living branch height (H b ), height at maximum canopy width (H c ), canopy width (C), age (Age), and height of trees (H) from Base_TAD with the growth model based on spatial structure [33] (GMBOSS) to obtain the current tree attribute data based on spatial structure (Base_TAD_SS).
We constructed a canopy curve using the under-living branch height, height at maximum canopy width, canopy width, and height of trees in Base_TAD_SS as the type value points of the B-spline curve, and the current canopy morphological data based on the spatial structure (Current_CMDD_SS) was obtained from the canopy shape curve.
The GMBOSS is expressed as Equations ( 8)- (10), where P v denotes the vertical spatial structure parameter, and P h denotes the horizontal spatial structure parameters.
Current_CMDD, Current_CMDD_CCLM, and Current_CMDD_SS are all canopy morphological data, which can be represented by a set (i, d, h, w).The distance is calculated with Equation (10), where x i and y i are the canopy widths at the corresponding locations, and n denotes the number of canopy description data.
We calculated the Euclidean distance D_CCLM between Current_CMDD and Cur-rent_CMDD_CCLM and the Euclidean distance D_SS between Current_CMDD and Cur-rent_CMDD_SS, and we compared the magnitudes of D_CCLM and D_SS to determine the degree of superiority of the two fitting effects.The contrast validation model is shown in Figure 6.

Calculation of PAR for Fir Tree Canopies Based on Ray-Tracing Algorithm
The annual average PAR in You County, Hunan Province, was 119.13 and 121.62 W/m 2 in 2015 and 2016, respectively.The daily variation in PAR showed a single-peak curve, which gradually increased from morning and reached the daily maximum around noon, and then began to gradually decrease.We combined the annual average PAR and the daily variation curve of the PAR, and we estimated the annual average values of the seven moments.The PAR intensity at each measurement point was obtained with the ray-tracing algorithm to calculate the summation to represent the PAR intensity at that point.A measurement point was taken at one-meter intervals for each fir tree in the four directions of the crown shape, and the PAR intensity (LDD) of the measurement point was calculated for a total of two years.The specific results are shown in Table 2.The maximum value of PAR at each point roughly positively correlated with the annual average PAR of the year.Due to the presence of forest stand closure, some points that receive little sunlight are always present.The average value of PAR at each point showed a decreasing trend, with little difference in the annual average PAR between the two years.Perhaps with the growth of the forest stand, the forest stand gradually becomes closed and light is blocked at more points in the forest stand, causing a decrease in the mean value.We used the 2015 survey data as Base_CMDD and the 2016 survey data as Cur-rent_CMDD; at present, the LDD is LDD15.All these data were used as the basis to obtain Current_CMDD_N with the BDGM and HDCM combined.The specific results are shown in Table 3: the maximum, minimum, mean, and standard deviation of Current_CMDD were 3.60, 0.20, 1.97, and 0.84, respectively.The maximum, minimum, mean, and standard deviation of Current_CMDD_N were 3.72, 0.15, 1.89, and 0.79, respectively.Current_CMDD describes the fir canopy morphological data from field surveys in 2016.We used Current_CMDD_N for the fir canopy morphological data from BDGM and HDCM simulations in 2016, which could indirectly represent the age factor.Each measurement point had a corresponding height, but the height was not used as a model parameter because the canopy width of the same height could be different for different fir trees.We used the multiple stepwise regression method, and Current_CMDD, Current_CMDD_N, and LDD were analyzed in SPSS.The results are shown in Table 4.With the introduction of Current_CMDD_N and LDD as the independent variables to the model, the fit of the model increased to some extent, with the coefficient of determination (R 2 ) increasing from 0.813 to 0.829, the corrected coefficient of determination (Rc 2 ) increasing from 0.806 to 0.821, and the standard estimation error decreasing from 0.3823 to 0.3491.The results indicated that Current_ CMDD_N and LDD had a stronger influence on the actual canopy width at the measurement point.The model with a better fit was selected as the final model, with a larger Rc 2 and a smaller standard estimation error (SEE).
The specific model parameters are shown in Table 5.For the CCLM, the p-values for the intercept, Current_CMDD_N, and LDD were 0.34, 0.018, and 0.025, respectively.The p-values of the three variables entering the model were all less than 0.05, showing significant differences.The final analysis determined that the fir CCLM was Current_CMDD = −40.418+ 0.927 * Current_CMDD_N + 0.012 * LDD (12)

Model Testing
As described in the above sections, based on the data from field surveys in 2015 and 2016, Current_CMDD, Current_CMDD_N, and LDD were fitted to obtain the fir CCLM.We then used the data from the field surveys in 2016 and 2017 to validate the model.We used the 2016 survey data as Base_CMDD and the 2017 survey data as Current_CMDD, at which time the LDD was LDD16.The Current_CMDD_N could be obtained by using the BDGM and HDCM.The LDD and Current_CMDD_N were incorporated into the fir CCLM to obtain the Current_CMDD_CCLM.Using the GMBOSS, we obtained the Cur-rent_CMDD_SS.We performed the t-test on the results, and the results are shown in Table 6.The table shows that the results of both the Current_CMDD_CCLM and Current_CMDD_SS were larger than that of the Current_CMDD, but the result of Current_CMDD_SS was larger.The p-values were all greater than 0.05, indicating a lack of significant difference between Current_CMDD_SS and the measured value, with no significant difference between Current_CMDD_CCLM and the measured value.

Model Comparison
Current_CMDD_CCLM, Current_CMDD_SS, and Current_CMDD are canopy morphological descriptions, all of which are a series of sets of canopy widths at measurement points from corresponding fir trees, directions, and heights.The 3D spatial coordinates of the measurement points were obtained by the coordinates of the fir roots (X, Y, Z), the directions of the measurement points, and the canopy widths of the measurement points, which we represented in the 3D space (Figure 7).The Euclidean distance D_SS between Current_CMDD_SS and Current_CMDD was 23.944, and the D_CCLM between Current_CMDD_CCLM and Current_CMDD was 15.561.The smaller the distance, the more accurate the simulation.This showed that the simulation using the CCLM is better than that using the GMBOSS at the forest-stand scale.
Current_CMDD_CCLM and Current_CMDD was 15.561.The smaller the distance, the more accurate the simulation.This showed that the simulation using the CCLM is better than that using the GMBOSS at the forest-stand scale.We attempted to identify differences in the canopy width (w) between Current_CMDD_ CCLM and Current_CMDD (Figure 8a) and between Current_CMDD_SS and Current_CMDD (Figure 8b).The simulation effect shown in Figure 8a is better than that in Figure 8b, demonstrating that the simulation using the CCLM is better than that using the GMBOSS.
Forests 2023, 14, x FOR PEER REVIEW 14 of 20 We attempted to identify differences in the canopy width (w) between Cur-rent_CMDD_CCLM and Current_CMDD (Figure 8a) and between Current_CMDD_SS and Current_CMDD (Figure 8b).The simulation effect shown in Figure 8a is better than that in Figure 8b, demonstrating that the simulation using the CCLM is better than that using the GMBOSS.The height of the fir trees in the sample plot ranged from 9.1 to 18.5 m.We randomly selected one fir tree each with a height of 10 m, 12 m, 14 m, 16 m, and 18 m.The canopy curves of the fir trees were plotted according to Current_CMDD_CCLM, Cur-rent_CMDD_SS, and Current_CMDD in two quadrants (east-west and north-south), as shown in Figure 9.And according to the canopy curves, three-dimensional models of the fir trees is constructed (Figure 10).Because the tree heights in Current_CMDD_CCLM and Current_CMDD_SS were both obtained from BDGM and HDCM simulations, the tree heights obtained from the simulations using both methods were the same.The tree heights obtained from the simulations were slightly higher than the measured tree heights, except for the fir tree with a height of 10 m.
(a) (b) The height of the fir trees in the sample plot ranged from 9.1 to 18.5 m.We randomly selected one fir tree each with a height of 10 m, 12 m, 14 m, 16 m, and 18 m.The canopy curves of the fir trees were plotted according to Current_CMDD_CCLM, Current_CMDD_SS, and Current_CMDD in two quadrants (east-west and north-south), as shown in Figure 9.And according to the canopy curves, three-dimensional models of the fir trees is constructed (Figure 10).Because the tree heights in Current_CMDD_CCLM and Current_CMDD_SS were both obtained from BDGM and HDCM simulations, the tree heights obtained from the simulations using both methods were the same.The tree heights obtained from the simulations were slightly higher than the measured tree heights, except for the fir tree with a height of 10 m.
As for the under-living branch height, height at maximum canopy width, and canopy width, the values of under-living branch height of the 10 m fir tree in the west and north, the height at the maximum canopy width of the 12 m fir tree in the west, and the height at maximum canopy width of the 18 m fir tree in the east obtained using GMBOSS simulations were closer to the measured values.For the rest of the values, the results obtained using CCLM simulations were closer to the measured values.
For the fit of the overall canopy curve, the GMBOSS simulation used the underliving branch height, height at maximum canopy width, canopy width, and height as the type value points of the canopy morphology and lacked control over other parts of the canopy morphological, so the fit between the canopy morphology obtained from the CCLM simulation and the actual canopy morphology was better than that between the GMBOSS simulation and actual canopy morphology.Therefore, at the single-tree scale, the results obtained using CCLM simulations were better than those obtained using GMBOSS simulations in the overall canopy curve or canopy morphological descriptors such as under-living branch height, height at maximum canopy width, and canopy width.
rent_CMDD_SS, and Current_CMDD in two quadrants (east-west and north-south), as shown in Figure 9.And according to the canopy curves, three-dimensional models of the fir trees is constructed (Figure 10).Because the tree heights in Current_CMDD_CCLM and Current_CMDD_SS were both obtained from BDGM and HDCM simulations, the tree heights obtained from the simulations using both methods were the same.The tree heights obtained from the simulations were slightly higher than the measured tree heights, except for the fir tree with a height of 10 m.As for the under-living branch height, height at maximum canopy width, and canopy width, the values of under-living branch height of the 10 m fir tree in the west and north, the height at the maximum canopy width of the 12 m fir tree in the west, and the height at maximum canopy width of the 18 m fir tree in the east obtained using GMBOSS simulations were closer to the measured values.For the rest of the values, the results obtained using CCLM simulations were closer to the measured values.
For the fit of the overall canopy curve, the GMBOSS simulation used the under-living branch height, height at maximum canopy width, canopy width, and height as the type value points of the canopy morphology and lacked control over other parts of the canopy morphological, so the fit between the canopy morphology obtained from the CCLM simulation and the actual canopy morphology was better than that between the GMBOSS

Discussion
In the traditional forest growth simulation process based on the BDGM and HDCM, the changes in diameter at breast height and tree height are considered, whereas the changes in canopy morphology are not.However, the environment has an important influence on canopy morphology.Some researchers have considered the influence of spatial structure and visually simulated fir tree canopies based on spatial structure.One of the most important factors for plant growth is light, and the influence of spatial structure on plant canopy morphology largely occurs because spatial structure affects the light distribution in the plant canopy, which in turn affects the photosynthesis of the plant and thus the canopy morphology.In studies on the effect of light on plant canopy morphology, researchers have quantitatively examined the effect of light on canopy morphology at the scale of a single tree and qualitatively studied this at the scale of the forest stands by describing light intensity through words such as "stand gap", "full light", and "shade".However, quantitative research on the effect of light on the canopy morphology at the forest-stand scale is lacking.In this study, based on the BDGM and HDCM, we introduced a ray-tracing algorithm from the computer field to calculate the light distribution in each part of the canopy in a virtual stand and combined this with field survey data, which simplified the acquisition of canopy light distribution data in a forest stand.The light distribution data were then introduced into the BDGM and HDCM to construct the CCLM, which quantitatively evaluated the effect of light on the canopy morphology at the forest-stand scale.
The canopy width in Huangfengqiao Forestry Farm was measured at one-meter intervals in four quadrants to describe the canopy morphology; in addition, we measured the under-living branch height, height at maximum canopy width, and canopy width of the fir trees in the four quadrants in the sample plot.A 3D scene was constructed with the actual data, and a ray-tracing algorithm was used in the 3D scene to calculate the light distribution data in each part of the tree canopy.Using the stepwise regression method, the relationships between canopy morphological data and age (the canopy morphological data obtained from the BDGM and HDCM were used to represent age) and light distribution data were analyzed.A CCLM of the fir tree was constructed with a coefficient of determination (R 2 ) of 0.829.The model is expressed in Equation (12).Current_CMDD was the current canopy morphological data obtained from the field surveys, Current_CMDD_N was the canopy morphological data obtained from the BDGM and HDCM simulations without adding any environmental factors, and LDD was the light distribution data.The value of the intercept was −40.418 (p = 0.034).The value of the Current_CMDD_N was 0.927 (p = 0.018).The value of the LDD was 0.012 (p = 0.025).
The canopy data from the field surveys, CCLM simulation, and GMBOSS simulation were compared to evaluate the fit of the canopy simulation using the CCLM and GMBOSS.For the forest-stand scale, the Euclidean distance between the overall canopy data was calculated to evaluate the fit with actual canopy morphology.The results indicated a divergence of 23.944 in the Euclidean distance between the canopy data from the field surveys and that from the GMBOSS simulation.The Euclidean distance between the canopy data from the field surveys and that from the CCLM simulation was 15.561.This showed that the simulation results using the CCLM were 8.383 times more accurate than those obtained using the GMBOSS.The fitting effect of the CCLM was better than that of the GMBOSS at the forest-stand scale.
For the single-tree scale, we randomly selected one 10 m high fir tree, one 12 m high fir tree, one 14 m high fir tree, one 16 m high fir tree, and one 18 m high fir tree, according to the height range of fir trees.The canopy curves of the fir trees were separately plotted based on the canopy data from the field surveys, CCLM simulation, and GMBOSS simulation to compare and analyze the fitting effect of the canopy curves of single trees of different heights.Because the tree heights obtained from the CCLM and GMBOSS simulations were simulated by the BDGM and HDCM, the height obtained from the simulations using the two methods were the same.As for the under-living branch height, height at maximum canopy width, and canopy width, the values of under-living branch height of the 10 m high fir tree in the west and north, the height at the maximum canopy width of the 12 m high fir tree in the west, and the height at maximum canopy width of the 18 m high fir tree in the east obtained using GMBOSS simulations were closer to the measured values.For the rest of the values, the results obtained using the CCLM simulation were closer to the measured values.For the fit of the canopy curve of a single tree, the GMBOSS simulation used under-living branch height, height at maximum canopy width, canopy width, and height as the type value points of the canopy morphology, which lacked control over other parts of the canopy morphology, so the fit between the canopy morphological data from the CCLM simulation and the actual canopy morphological data was better than that between the canopy morphology from GMBOSS simulation and the actual canopy.Therefore, at the single-tree scale, the results obtained from the CCLM simulations were more accurate than those obtained from GMBOSS simulations in the overall canopy curve or canopy morphological descriptors such as under-living branch height, height at maximum canopy width, and canopy width.
This study is based on the BDGM and HDCM by introducing the factor of light to construct the CCLM, so the accuracy of the CCLM is increased based on the accuracy of the BDGM and HDCM.As such, the accuracy of the selected BDGM and HDCM strongly influenced on the accuracy of results.The light distribution data in this study were calculated by using the annual average PAR with the ray-tracing algorithm.Although the ray-tracing algorithm can accurately simulate realistic light situations, the study is not free from errors, which may also have had an impact on the experimental results.Other factors that may affect the canopy morphology include broken tips, pests, and diseases, thereby influencing the accuracy of the study results.Additionally, both tree growth and solar movement are continuously changing processes.In this study, the initial state of the stand was used as the basis, and the solar position at seven moments was selected to calculate the light distribution of the canopy.These factors deserve more attention in future research.

Conclusions
In this study, we considered the effect of light on the canopy shape and constructed a new coupled canopy-light model (CCLM) based on light distribution data and canopy morphological data, which was used to improve the visual polymorphism simulation of fir morphology.The model is expressed in Equation ( 12) with a coefficient of determination (R2) of 0.829.The value of the intercept was −40.418 (p = 0.034).The value of the Cur-rent_CMDD_N was 0.927 (p = 0.018).The value of the LDD was 0.012 (p = 0.025).For the forest-stand scale, the Euclidean distance between the canopy data from the field surveys and that from the CCLM simulation was 15.561.The CCLM was compared with the GM-BOSS proposed by Ma et al., and the Euclidean distance between the canopy data from the field surveys and that from the GMBOSS simulation was 23.944.The results indicated that the simulation results using the CCLM were 8.383 times smaller than those obtained using the GMBOSS.This means that the CCLM was 35.0% more accurate than the GMBOSS in the visual polymorphism simulation of fir morphology, and the fitting effect of the CCLM was better than that of the GMBOSS at the forest-stand scale.The computer graphics technology was closely combined with constructing tree canopy shape and improved visual polymorphism simulation of fir morphology, which will enable us to provide more accurate tree canopy models for forestry research.However, to make visual polymorphism simulation more accurate, we will consider the impact of other environmental factors in the future.

Figure 1 .
Figure 1.Description of study area: (a) Location of study area in Hunan Province.(b) Enlargement of inset in (a) showing details of the Huangfengqiao Forest and the sample area localization in Yu Country.(c) Example of the vegetation structure with forest canopy aspects.

Figure 2 .
Figure 2. Spatial distribution of fir trees in sample plot.The blue points are the coordinate points representing the relative position of the fir.

Figure 1 .
Figure 1.Description of study area: (a) Location of study area in Hunan Province.(b) Enlargement of inset in (a) showing details of the Huangfengqiao Forest and the sample area localization in Yu Country.(c) Example of the vegetation structure with forest canopy aspects.

ForestsFigure 1 .
Figure 1.Description of study area: (a) Location of study area in Hunan Province.(b) Enlargement of inset in (a) showing details of the Huangfengqiao Forest and the sample area localization in Yu Country.(c) Example of the vegetation structure with forest canopy aspects.

Figure 2 .
Figure 2. Spatial distribution of fir trees in sample plot.The blue points are the coordinate points representing the relative position of the fir.

Figure 2 .
Figure 2. Spatial distribution of fir trees in sample plot.The blue points are the coordinate points representing the relative position of the fir.

Figure 3 .
Figure 3. Description of canopy shape and evaluation of fit.The x and y represent the distance from the canopy edge to the trunk; the arrows with two directions are used to identify this distance.

Figure 3 .
Figure 3. Description of canopy shape and evaluation of fit.The x and y represent the distance from the canopy edge to the trunk; the arrows with two directions are used to identify this distance.

Figure 4 .
Figure 4. Schematic diagram of ray-tracing algorithm: (a) light starts from light source and eventually hits a camera; (b) light shines out of the camera and backtracks.

Figure 4 .
Figure 4. Schematic diagram of ray-tracing algorithm: (a) light starts from light source and eventually hits a camera; (b) light shines out of the camera and backtracks.

Table 1 .
Basic information of fir trees in sample site.

Table 1 .
Basic information of fir trees in sample site.

Table 2 .
Photosynthetically active radiation (PAR) at measurement points in 2015 and 2016.

Table 3 .
Canopy morphological data from field surveys and from BDGM and HDCM simulations in 2016.

Table 5 .
Statistics of the CCLM.

Table 6 .
Results of t-test of CCLM and GMBOSS simulations.