1.

Introduction

Plant’s growth and development can be negatively affected by biotic or abiotic stress factors.¹ Biotic stress is caused by other organisms, such as insects, pathogens, or other plants. Abiotic stress refers to physical and chemical factors in growth environment, such as pollution, drought, and nutrient availability. Remotely sensing vegetation condition and health by hyperspectral imaging (imaging spectroscopy) has been widely studied in recent decades.² High spectral resolution of hyperspectral images enables finding spectral bands that are most distinctly affected by a certain stressor.² In forestry, hyperspectral remote sensing has been used to detect pests and pathogens,^3–5 ozone damage,⁶ effects of pollution,⁷ and identifying tree species,⁸ for example.

Analysis of remotely sensed spectral images relies largely on machine learning (ML) algorithms, although statistical tests have also been widely applied.² As training of ML algorithms requires large amounts of training data cumbersome to gather from real-world measurements, it is often generated with simulators. In practice, the simulators, or leaf optical properties models, mimic spectral response of plants with known biophysical and biochemical properties. Predicting the presence and type of stressors from a real hyperspectral image is done by finding spectral features that match the simulation, i.e., the simulation problem is inverted.

A rough taxonomy of the leaf optical properties models can be established by examining the scale that they operate on: cell scale models depict the intricate internal structure of a plant leaf,^9–11 leaf scale models see the leaves as a whole,^12,13 and canopy scale models simulate the reflectance of groups of plants.^14–17 Canopy scale simulators are often statistical extensions of leaf models with some additional elements, such as soil reflectance. Canopy scale models can simulate remotely sensed hyperspectral images and thus, they are essential in training ML algorithms.

Most existing canopy scale simulators are developed to simulate high altitude imaging where an imager is mounted to a satellite or an airplane. The first satellite with spectral imaging equipment was the Landsat 1 mission (originally called Earth Resources Technology Satellite, ERTS) launched in 1972 that was equipped with four-band multispectral imager.¹⁸ The first airplane-based hyperspectral imaging sensor capable of measuring continuous spectra with 10 nm bandwidth in the 400 to 2500 nm range commonly used in vegetation spectroscopy was the airborne visible/infrared imaging spectrometer (AVIRIS).¹⁹ AVIRIS produced its first hyperspectral images in 1987.¹⁹ With growing popularity and decreasing prices of unmanned aerial vehicles (commonly called drones) and hyperspectral imagers, there is need for leaf optical simulators aimed at low altitude imaging. For a comprehensive review of drone-based imaging, see Ref. 20.

One of the most noticeable differences between high and low altitude imaging is in the ground pixel size (spatial resolution), i.e., how many meters in ground a single pixel covers. The ground pixel size of airplane-based imagers range from 0.5 to 10 m and is even greater in space-borne imagers.² Utilizing statistical approximations for canopy structure, such as the commonly used average leaf inclination angle (LIA), is well-founded in these applications. However, in low altitude drone-based imaging, ground pixel size may be centimeters or even millimeters depending on the sensor and altitude. As such, the concept of average LIA is not applicable since a pixel can be filled by single leaf and the inclination angle of that leaf is what it is. Simulating low altitude imaging demands intricate description of plant geometry.

One of the most recent attempts to address the need for a low altitude vegetation simulator is project HyperBlend.²¹ HyperBlend’s working principle is shown in Fig. 1. The simulator requires reflectance and transmittance spectra of a leaf as an input; they are called target spectra. The target spectra can originate from real-world measurements, or they can be generated by a third party leaf simulator. In the following conversion step, HyperBlend creates a renderable three-dimensionl (3D)-model of a leaf that respects spectral reflectance and transmittance properties defined by the target spectra. Copies of the leaves are then placed on 3D-models of trees, which in turn are placed on 3D-model of a forest floor to form a forest canopy model. In Fig. 2, an example of a tree generated by this process is shown. Finally, a virtual camera renders an image cube of the forest accompanied by pixel labels that have information on which object was present in each pixel.

Fig. 1

Working principle of HyperBlend canopy model. Steps considered in this paper are highlighted in red.

Fig. 2

An example of a tree generated with HyperBlend. (a) Pure geometry without shading. (b) Tree model shaded with arbitrary RGB colors.

Rendering is a computer graphics term that means the process of producing a synthetic image of a virtual scene.²² HyperBlend utilizes 3D-modeling and rendering software Blender²³ and especially its path tracing rendering engine Cycles. Path tracing is a rendering technique that numerically approximates Kajiya’s rendering equation.²⁴ Kajiya’s path tracing only accounts for light interactions on object’s surface.²⁴ Later efforts of Lafortune extended the algorithm to include interactions inside an object as well.²⁵ HyperBlend takes advantage of volumetric path tracing by depicting a leaf as translucent media that can scatter and absorb light.

In this study, we only consider HyperBlend’s leaf scale operations, i.e., we exclude tree and forest geometry generation and construction of the image cube. HyperBlend’s leaf model depicts a leaf as a flat box whose volume is filled with absorbing and scattering particles. Desired reflectance and transmittance properties can be achieved by adjusting absorbing particle density, scattering particle density, forward versus backwards scattering tendency, and mixing factor of the two particle types. For more detailed explanation of the working principle of HyperBlend’s leaf model, see Ref. 21. HyperBlend is written in Python 3.7²⁶ and the source code is available in a GitHub repository at: https://github.com/silmae/hyperblend. The version that this paper refers to is tagged with version number 0.2.0.

While the conversion step (highlighted in red in Fig. 1) in the original version of HyperBlend produced highly accurate results with root mean squared error (RMSE) of $\sim 1\%$,²¹ in this paper, we wish to address some of its limitations: simulation speed, parameterization, and generalizability.

The conversion step of the original version was very slow taking more than an hour to convert a pair of reflectance and transmittance spectra into renderable leaf material parameters on a consumer-grade computer.²¹ We present two new methods for the conversion that are computationally more efficient: a surface fitting method in Sec. 2.2 and a neural network (NN) method in Sec. 2.3. Training data generation for the new methods is explained in Sec. 2.1.

HyperBlend relies on an external source of target spectra (top left corner in Fig. 1). Target spectra should include associated biophysical and biochemical parameters, such as water content and chlorophyll content, which explain the shape of the spectra. Although real-world datasets where measured spectra are accompanied by biophysical and biochemical data exist, it is often more convenient to produce such data with a leaf simulator. As HyperBlend’s conversion step is agnostic to the source of such data, both methods can be used. In Sec. 2.4, we include an implementation of the well-known PROSPECT^13,27–31 leaf simulator into HyperBlend. The role of the PROSPECT model in HyperBlend is the same as, for example, in PROSAIL^17,29 where PROSPECT acts as an underlying leaf model that is combined with canopy model SAIL.¹⁴ From the perspective of the user of HyperBlend, this also hides the unwieldy rendering parameters under much simpler PROSPECT parameterization.

Lastly, as the validation of the original model in Ref. 21 utilized leaf samples of only a single species, HyperBlend’s generalizability was not verified. In this study, we run PROSPECT with randomized parameters to generate a wide variety of target spectra to show that HyperBlend can accurately reproduce leaf spectra also in a general case.

We conduct two simulation experiments in Sec. 3. The first one is a repeat experiment of Ref. 21, which compares simulation speed and accuracy of the presented methods against the original simulation method. The second experiment tests how well HyperBlend can reproduce arbitrary leaf spectra. Results of the experiments are presented in Sec. 4. Future developement and limitations are discussed in Sec. 5 and Sec. 6 concludes the paper.

2.

Materials and Methods

The leaf model in HyperBlend works by adjusting four leaf matter parameters (absorbing particle density ${\rho }_a$, scattering particle density ${\rho }_s$, scattering anisotropy $\alpha$, and mixing factor of the two particle types $\beta$) until target reflection and transmission values are achieved. The target can be obtained from real world measurements or from another simulator. As Blender is not capable of working with spectral quantities, HyperBlend simulates each wavelength band separately before compiling them into a spectrum. Thus, we are mainly interested in which values of ${\rho }_a$, ${\rho }_s$, $\alpha$, $\beta$ can produce a certain reflectance-transmittance pair (RT pair).

Letting $R$ be reflectance in closed interval [0, 1] and $T$ transmittance in the same interval, we can denote the four leaf material parameters by four respective functions

In the rest of the paper, we will use three different methods in running HyperBlend’s leaf model. The first is the original optimization-based method presented in Ref. 21 which we will call the “Original” method. Two new methods are presented below: “Surface” fitting method and NN method. These short names will be used for brevity throughout the rest of the paper.

All computations were done with a consumer-grade machine with Intel i9-9900K 8-core processor. Hereafter, all the plots in this paper are plotted with Matplotlib.³²

2.1.

Training Data Generation

In order to find functions $f_i$ in Eq. (1), we generated evenly spaced RT pairs with following rules:

Eq. (2)

$$\begin{gathered} |r-t|\le 0.25 \\ r+t\le 1 \\ r,t\ge 0. \end{gathered}$$

The first rule ensures that reflectance and transmittance are always close to each other as is usually the case in plant leaves. The threshold value of 0.25 was selected because the accuracy of the Original optimization method begins to decrease when reflectance and transmittance values are too dissimilar. The second rule ensures that no more photons can be reflected and transmitted than what is incident on the leaf, i.e., bioluminescence is not allowed. The last rule ensures that reflectance or transmittance cannot be negative.

After generating the RT pairs, Original HyperBlend leaf model was run for each of them to produce leaf material parameters ${\rho }_a$, ${\rho }_s$, $\alpha$, $\beta$, which, in turn, produce simulated reflectance $\widehat{r}$ and simulated transmittance $\widehat{t}$. To avoid grossly bad data points, a point was removed from the training set if either $|r-\widehat{r}|\ge 0.01$ or $|t-\widehat{t}|\ge 0.01$.

From initial point generation runs, we noticed that points where $t>3.8r+0.02$ or $t<0.5r-0.035$ were almost always exceeding the pruning limit. We added those areas as additional conditions where training points were not generated to avoid unnecessary computing of points that would be pruned later. Figure 3 shows generated training data, the additional conditions, and pruned points. One can notice from the figure that accuracy of the simulation decreases toward the edges of the area. The final number of accepted training points was 44,284. The training data set is available at https://osf.io/trhf8/.

Fig. 3

Generated training data points. Black lines show additional restrictions. Points shown in red exceed the 0.01 error threshold and were pruned from the final training set. Only points shown in blue were used in final training. Axes $R$ and $T$ stand for reflectance and transmittance, respectively.

2.2.

Surface Fitting

Now that we have simulated leaf material parameters for each RT pair, we can plot them into $R$, $T$, $z$-space, where $z$-axis shows the value of each leaf material parameter for given $r$ and $t$ (see Fig. 4). We assumed that the leaf material parameters behave smoothly in that space so we can fit an analytical function to approximate each material parameter.

Fig. 4

Leaf material parameter surfaces shown as training data (light gray), surface fitting (orange), and trained NN (green). The differences between surface fitting and NN are easiest to see at the bottom right corner of the plots. For all parameters, the NN fits the training data better than the surface fit. (a) Absorption density ${\rho }_a$, (b) scattering density ${\rho }_s$, (c) scattering anisotropy $\alpha$, and (d) mix factor $\beta$.

The first suggested method to find functions $f_i$ in Eq. (1) is to fit a surface into training data. As shown in Fig. 4, the shape of the surfaces varies. We used exponential, logarithmic, and polynomial functions depending on the shape of the surface:

Eq. (3)

$$\begin{gathered} {\rho }_a=f_0(r,t)=a_0{\mathrm{e}}^{b_{0}r}+c_0{\mathrm{e}}^{d_{0}t} \\ {\rho }_s=f_1(r,t)=a_1\log (b_{1}r)+c_1\log (d_{1}t) \\ \alpha =f_2(r,t)=a_2{r}^{b_2}+c_2{t}^{d_2} \\ \beta =f_3(r,t)=a_3{\mathrm{e}}^{b_{3}r}+c_3{\mathrm{e}}^{d_{3}t}, \end{gathered}$$

where $a_i,b_i,c_i,d_i$ are constants to fit each function into the training data. In this case, function fitting is a non-linear least squares problem that minimizes the squared $L1$ distance of the prediction function $f_i(r,t)$ to the training data $z_i(r,t)$ at the same point

Eq. (4)

$$\underset{a_i,b_i,c_i,d_i}{\min }\sum _{k=1}^n{(z_i(r_k,t_k)-f_i(r_k,t_k))}^2,$$

where $n$ is the number of training points, subscript $i$ refers each of the four leaf material parameters. This kind of fitting problem can readily be solved with Levenberg–Marquardt algorithm.^33,34 For the code implementation, we used SciPy’s³⁵ optimize.curve_fit method with the default options.

The surfaces align fairly well with most of the training data but problems do occur especially when both reflectance and transmittance are close to zero. This behavior is especially apparent in the plot of anisotropy parameter $\alpha$ in Fig. 4(c) where the surface twists into almost vertical orientation near the origin. The training data also have ripples in the direction of $r=t$ as shown in Fig. 4. These features appear too regular to be random noise, suggesting that they may contain some relevant information.

3.

Experiments

To assess the performance of the new Surface and NN simulation methods, we conducted two experiments. The first experiment replicates the one in the original paper.²¹ The second experiment tests HyperBlend’s capability to generalize plant leaf spectra with varying characteristics. Both experiments were run with 5 nm spectral resolution.

Experiment 1 is essentially the same experiment that was conducted in the original HyperBlend paper.²¹ The only difference is that now, in addition to the Original method, we also use presented Surface and NN methods. We use the exact same data set of 28 spectral measurements of silver birch (Betula pendula) leaves collected by Ref. 42. For more details on the original experiment, see Ref. 21. The main reason for this experiment is to compare the simulation time and accuracy of the old method to the new methods presented above.

Experiment 2 tests generalizability of HyperBlend. HyperBlend has previously been tested only against the data set used in experiment 1 that consists of a samples of single tree species, so HyperBlend’s generalizability has not previously been demonstrated. This experiment also shows that the parameter conversion from PROSPECT parameters is valid. We assessed the accuracy of conversion by generating 100 random leaf spectra with PROSPECT. For random generation of the spectra, we used the standard PROSPECT leaf of Ref. 43 (p. 275) as a basis, and let each parameter vary independently by drawing from normal distribution centered at the standard value. Drawn values were clipped to range suggested by Ref. 44 whose PROSPECT implementation we adapted to our code base. Mean and variance for normal distributions along with clipping ranges for each PROSPECT parameter are shown in Table 1. Experiment 2 was not run with Original method due to the excessive computation time this would have required.

Table 1

PROSPECT spectra were generated by drawing from normal distribution using mean and variance tabulated here. If necessary, drawn values were clipped to the shown range.

4.

Results

Mean errors and simulation times for both experiments are shown in Table 2. Mean errors are shown over all samples and over the whole spectrum. Simulation time is the mean over all samples.

Table 2

Comparison of simulation speed and errors between Original method, Surface method, and NN method. Simulation speed is reported in seconds per sample and errors for reflectance (R) and transmittance (T) as RMSE. Experiment 2 was not run with Original method. Simulation times for Surface and NN are roughly 200 times smaller than that of the Original. While both sacrifice some accuracy, NN is consistently more accurate than Surface.

In the first experiment with birch leaf data, both new methods are roughly 200 times faster than the Original with 11 to 12 s per sample. For the whole 28 sample set in experiment 1, this means reducing simulation time from 16 h to less than 6 min. Simulation time of presented methods stay consistent in experiment 2, which shows that simulation speed is constant regardless of the shape of the simulated spectra. This is not true with the Original method, as the optimization process takes more time if the starting guess is not close to the optimal solution.

The cost of faster simulation is increase in simulation error. RMSE of the Surface method in experiment 1 is roughly 10 times larger than the Original method. The error in NN is roughly doubled. We consider doubling the error for 200 times speedup a good trade. The errors of experiment 2 show that the accuracy of NN method does not suffer even when more difficult data are used. The Surface method, on the other hand, has difficulties especially in simulating transmittance.

More detailed spectrum-wise results of the experiment 1 are shown in Fig. 5. Even if the mean error 0.01 of the Surface method does not seem much, it is clearly visible in several areas of mean spectrum in Fig. 5(b). Large difficulties can be seen in short wavelengths ($<500\mathrm{nm}$) and in near infrared (NIR) plateau (700 to 1300 nm) on both reflectance and transmittance and in long wavelengths ($>2100\mathrm{nm}$) transmittance. The NN method in Fig. 5(c) does not show similar difficulties.

Fig. 5

Mean spectra (top row) and mean RMSE (bottom row) of the 28 silver birch samples of experiment 1. In the top row plots, target mean spectra plotted as black dotted line, target variance as gray dashed line, simulated mean as colored solid line, and simulated variance as colored shadow. Simulation error of Surface method is one order of magnitude greater than that of the other two methods. (a) Original, (b) Suface, (c) NN, (d) Original, (e) Surface, and (f) NN

Similar plots for experiment 2 are shown in Fig. 6. The accuracy of the NN method is similar to experiment 1, even though some over-estimation in transmittance can be seen in the NIR plateau in Fig. 6(b). The Surface method, on the other hand, has great difficulties in estimation of transmittance throughout the spectrum shown in Fig. 6(a). Estimate of the reflectance is nearly perfect with NN and on par with experiment 1 with Surface.

Fig. 6

Mean spectra (top row) and mean error (bottom row) of the of 100 PROSPECT-generated random leaf spectra samples of experiment 2. Target (i.e., PROSPECT-generated) mean spectra plotted as black dotted line, target variance as gray dashed line, simulated mean as colored solid line, and simulated variance as colored shadow. Simulation errors of both methods are slightly higher than in experiment 1. Again, the error of the Surface method is one order of magnitude higher than that of the NN method. (a) Surface, (b) NN, (c) Surface, and (d) NN.

To get more insight on why the errors are greater in experiment 2 than in experiment 1, we examined the solitary samples. Simulation results for sample number 59 are shown in Fig. 7. Inspecting the target spectra shown in black, one notices that the difference of reflectance and transmittance is quite large. For example, on wavelength 1000 nm in NIR plateau, the reflectance is 0.26 and transmittance is 0.60. The difference 0.34 is much greater than 0.25 that was allowed in our training data. Lack of training data depicting this kind of situation explains why Surface and NN perform poorly with this and seven other samples with similar difference in reflectance and transmittance.

Fig. 7

Difficult PROSPECT-generated sample number 59 simulated with all three methods. Black points show the target spectra and colored ones the simulation results. This sample shows the instability of our model if difference between reflectance and transmittance is too high. The Original method is unstable in the NIR plateau and systematically overestimates reflectance on long wavelengths. Training data of Surface and NN did not contain a sample with a similar situation, so they are merely guessing the spectra. (a) Original, (b) Surface, and (c) NN.

The similarity condition 0.25 was set in Eq. (2) because most of the erroneous points in generated training data were concentrated on the areas of greatest dissimilarity between reflectance and transmittance. For further verification that the similarity condition is indeed crucial, we simulated sample 59 with the Original method. The result of that simulation in Fig. 7(a) clearly shows that the Original method is unstable in the NIR plateau where some wavelengths are very poorly simulated while others are at least approximately correct. Based on that behavior, it seems that the optimal solution does indeed exist, but it has not been reached in all cases. On the other hand, looking at long wavelengths, correct reflectance is never reached so the solution might not exist at all. Both Surface in Fig. 7(b) and NN in Fig. 7(c) can find approximately correct reflectance spectra but fail spectacularly with transmittance. The other seven poorly simulated spectra (not shown) exhibit the same behavior.

We assume that the reason why simulation of dissimilar reflectance and transmittance fails has to do with scattering anisotropy parameter of cycle’s volume scattering shader. We have noticed that it can behave discontinuously when too small or too great values are used. Scattering anisotropy is mainly responsible for providing required asymmetry between reflectance and transmittance. For example, increasing scattering anisotropy favors forward scattering allowing more photons to pass the virtual leaf, which increases transmittance and decreases reflectance.

5.

Discussion

With the two new methods presented, HyperBlend can now convert spectra to renderable leaf parameters with three different methods: Original optimization method, Surface fitting method, and NN method. If maximal accuracy is needed, we recommend using Original and for vastly faster evaluation when slight inaccuracy can be allowed, we recommend using NN. We do not recommend using Surface at all because while it is as fast as NN, its accuracy is significantly worse.

Regardless of the method, the most noticeable error is made when reflectance and transmittance are too dissimilar. This should not affect modeling real leaves too much as they tend to have fairly similar reflectance and transmittance spectra. As Ref. 43 points out, drawing PROSPECT parameters as independent variables does not represent natural interconnections of biophysical and biochemical parameters, which may result in spectra that does not occur in real leaves.

Regarding the failure to simulate dissimilar reflectance and transmittance with the presented methods, one solution could be to generate more dissimilar training data. As the Original method often makes gross errors in these situations, a lot of points would have to be generated to get sufficient amount of usable data. Also, alleviating the 1% error threshold for acceptance should be considered, as even the best of dissimilar points are not very accurate. On the other hand, using less accurate training data would probably worsen the performance of the models in general. Generating a lot of training data takes a lot of time. Generating the 44,284 training points used in this study took roughly 72 h with a consumer-grade machine with Intel i9-9900K 8-core processor.

The leaf model could still be improved by adding surface reflections, which would produce more realistic leaves and allow modeling some special cases such as waxy leaves. Using bidirectional reflectance distribution function and bidirectional scattering distribution function for the leaf surface would be possible in Blender when those functions are known from real-world measurements. Another improvement would be to model needle-like leaves, which would make modeling coniferous and mixed forests possible.

6.

Conclusions

With the developments presented in this study, the conversion from target spectra to virtual renderable leaf material can now be done much faster in HyperBlend. While the conversion with the new methods is not as accurate as the original optimization method, the 200-fold increase in simulation speed makes HyperBlend more usable especially in quick drafting.

The new methods also produce more stable spectra. The original method is based on optimization, which has no guarantees to converge in a specific solution and may result in discrepancy in adjacent wavebands. As the functions used in the surface fitting method are smooth by definition, the resulting spectra are guaranteed to be smooth. The NN method does not have this kind of inbuilt guarantee, but the network trained in this study appears to be smooth enough to alleviate any concerns of jagged spectra.

HyperBlend is not a leaf optical simulator in the traditional sense that it could take some biophysical and biochemical parameters, such as chlorophyll content, as an input, and produce leaf reflectance and transmittance spectra as an output. However, by incorporating PROSPECT leaf optical model as a source of target spectra, from the user’s perspective, the functionality is as if they were running PROSPECT simulator in a 3D virtual world. Because the underlying leaf simulator is not tightly connected to 3D-graphics, it can be easily replaced with another simulator or measured data if needed. This gives a high degree of flexibility in adjusting the configuration to fit different purposes.

Future studies will cover generation of tree and forest geometry as well as construction of hyperspectral image cubes and associated metadata needed for inversion of the canopy model. This and the previous study on HyperBlend have focused on necessary prerequisites for creating simulated hyperspectral images of forests with Blender. The following studies will highlight the advantages of our image generation scheme: procedural geometry, largely automated workflow, fast rendering of complex scenes, and vast heterogeneous forest canopies complete with understory.