Efficient Storage and Importance Sampling for Fluorescent Reflectance

We propose a technique for efficient storage and importance sampling of fluorescent spectral data. Fluorescence is fully described by a re‐radiation matrix, which for a given input wavelength indicates how much energy is re‐emitted at other wavelengths. However, such representation has a considerable memory footprint. To significantly reduce memory requirements, we propose the use of Gaussian mixture models for the representation of re‐radiation matrices. Instead of the full‐resolution matrix, we work with a set of Gaussian parameters that also allow direct importance sampling. Furthermore, if accuracy is of concern, a re‐radiation matrix can be used jointly with efficient importance sampling provided by the Gaussian mixture. In this paper, we present our pipeline for efficient storage of bispectral data and provide its extensive evaluation on a large set of bispectral measurements. We show that our method is robust and colour accurate even with its comparably minor memory requirements and that it can be seamlessly integrated into a standard Monte Carlo path tracer.


Introduction
Traditionally, rendering pipelines mimic the human visual system and rely on RGB or tristimulus values to compute light transport. This approach, however, is insufficient if high colour accuracy or reproduction of spectral effects is required. Resolving the wavelength domain with additional bands, and thus performing spectral rendering, allows handling of such cases at the expense of increased computational demands. However, due to the increase in computation power and algorithmic efficiency [WND*14], spectral rendering is becoming increasingly popular even in production [PBC*18]. In this paper, we focus on the reproduction difficulties of one particular spectral effect: fluorescence.
Fluorescence is a phenomenon often found in natural elements and chemical compounds. It has a significant impact on how we perceive an object's colour in terms of hue, saturation and luminance because it transfers a part of non-visible or near-visible light into a different part of the visible range. Therefore, the inclusion of fluorescence can enhance material brightness, making the colour appear more vivid compared to a pure reflection. While the importance of this effect for appearance modelling in computer graphics has been known for years, it has not been used much in practice, as adding it to graphics workflows-both on the side of user interfaces, as well as within rendering software proper-is a demanding problem.
An issue preventing a more widespread adoption of fluorescence is the significant memory footprint of its representation. To fully describe fluorescent spectra, one typically uses re-radiation matrices [Don54], also called Donaldson matrices, representing reemission distributions for a given set of wavelengths. Although explicitly storing these data is manageable for a few fluorescent elements, this does not scale for fine-grained material definitions such as fluorescent textures.
In recent years, there has been an increasing interest in fluorescence in the research community. For example, Jung et al. [JWH*19] used fluorescence to enlarge the colour gamut during spectral uplifting in a fashion analogous to the use of optical brighteners in real-world materials. However, due to the lack of  (left), fluorescent surfaces also re-emit a portion of the absorbed energy from incident light as a light at additional wavelengths, leading to the colourful appearance of other objects in the scene. To represent fluorescent materials in a renderer, we typically use re-radiation matrices (centre), which have a significant memory overhead. Instead, we propose a more efficient representation of these matrices using Gaussian mixture models (right -eight Gaussians). Although compact, this representation also provides compelling results even with such challenging illumination.
better alternatives, their method had to store a full re-radiation matrix for every texel, which in practice limited its usage to low-resolution input images.
We have proposed a method for efficient storage and importance sampling of re-radiation matrices using a Gaussian mixture model ( Figure 1) in Hua et al. [HFW21]. This preliminary work was evaluated on three measured fluorescent materials. In this revision, we present an extended evaluation using an additional database of measured materials [GF00]. Further, we compare the previously proposed weighted Expectation-Maximization (EM) method with the weighted Bayesian approach of Gaussian mixture fitting. Finally, we provide additional evaluation metrics to assess the fidelity of the fitted re-radiation matrices.

Bispectral rendering
Although not widely used, the concept of fluorescence in computer graphics is not new. Glassner [Gla95] introduced fluorescence and phosphorescence to the computer graphics community and proposed a method to support those effects in a Whitted ray tracer. Wilkie et al. [WTP01] proposed a path tracing system capable of rendering both fluorescence and polarization effects, but they relied on a defined number of spectral bands, the wavelength not being part of the (MC) integration. Mojzík et al. [MFW18] made the wavelength domain part of the integrand and adapted Hero Wavelength Spectral Sampling (HWSS) [WND*14] to their method capable of handling fluorescent media in a Monte Carlo (MC) path tracer.

Colour gamut enlargement
One of the challenges when using a spectral renderer is the requirement of spectrally defined assets. Assets are commonly defined in a tristimulus space (e.g., RGB, XYZ) and cannot be used directly as the input of a spectral renderer. To allow the usage of RGB assets in a spectral renderer, the original RGB data needs to be 'uplifted' to a spectral form. A wide range of work addresses nonfluorescent spectral uplifting, assuming smooth spectra [Smi99, [JH19], or data-driven from measured spectra [OYH18,TWF21]. Recently, Jung et al. [JWH*19] proposed an uplifting pipeline that additionally uses fluorescent spectra to enhance the colour gamut. A significant drawback of this bispectral uplifting method is its memory requirement. In the worst case, when dealing with texture uplifting, each pixel needs to be stored as a full re-radiation matrix. However, this is not an intrinsic drawback of the uplifting method, but an indicator that the representation of fluorescent data in rendering technology required improvement.

Bispectral measurements
Bispectral reflectance measurements are not as widespread as BRDF measurements because their acquisition is a lengthy process requiring expensive equipment. Additionally, with the limited number of renderers supporting fluorescence, there is little incentive to take bispectral material measurements for computer graphics applications. However, few exceptions are noteworthy. Gonzalez and Fairchild [GF00] offer a database of re-radiation matrices of numerous materials, mainly papers and inks. They used a Labsphere BFC-450 bispectral fluorescence colorimeter for the measurements. We use this database to consolidate our earlier published evaluation on the three re-radiation matrices from Labsphere distributed in ART [ART18] measured with the same device.

Bispectral material models
Wilkie et al. [WWLP06] proposed one of the early bispectral models adapted to current rendering techniques. This model uses a layered BRDF with a diffuse fluorescent component. Hullin et al. [HHA*10] proposed an efficient acquisition setup to capture bispectral bidirectional reflectance and re-radiation distribution functions (BRRDF). They guided their acquisition by a principal component analysis to lower the acquisition time. Based on this work, Jung et al. [JHMD18] derived a new BRRDF. They exploited Kasha's rule to represent their distribution with three 1D distributions (absorption and emission spectra and non-fluorescent reflectance) and two ratios (emitted to absorbed energy and reradiation to non-fluorescent reflectance).

Gaussian representations in rendering
Gaussian distributions are often used in rendering practice, as their sufficiently large combination can approximate distributions of any shape. Jakob et al. [JRJ11] represent radiance using a Gaussian mixture and applied it to the rendering of volumetric media.

Background
In this section, we review the core concepts used in the remainder of this article. We explain the fluorescent effect, its properties and the typical representation of fluorescent data in a bispectral renderer. We also review the basic notation of Gaussian mixture models (GMMs), which we employ to represent re-radiation matrices.

Fluorescence
Fluorescence is an effect where a fluorophore absorbs energy from incident photons, causing excitation of its molecules which, during their return to the ground state, may re-emit part of this energy as photons with lower energy. In contrast to phosphorescence, where emission can transpire over a long time period, fluorescence occurs in a few nanoseconds (10 −9 to 10 −7 s) and, as such, is considered to be instantaneous in the context of computer graphics.
Fluorescence produces noticeable visual effects and is commonly used to increase the brightness and saturation of pigments and dyes. This works by shifting light from the barely visible ultra-violet region to longer wavelengths, where the human eye is more sensitive.
A fluorophore can be characterized by its absorption and reemission spectrum (Figure 2a). The first defines which incident wavelengths are absorbed and lead to re-emission events, with magnitudes representing the total intensity of the re-emissions. The second describes the amount of re-emission across all incident wavelengths.
According to Kasha's rule, the re-emission spectrum shape can be considered independent of the specific incident wavelength that caused the excitation (re-emission intensity scales). The difference between the spectral positions of the band maxima of absorption and re-emission arising from the same electronic transition is called a Stokes shift [MW97]. In practice, this shift usually corresponds to Figure 2: (a) If Kasha's rule holds, the shape of a re-emission spectrum is the same regardless of the incident excitation wavelength, and absorption with re-emission spectrum form a sufficient representation of a fluorophore. (b) Otherwise, a re-radiation matrix providing precise information about the re-emission spectrum given an incident wavelength has to be used at the cost of its significant memory footprint.
a difference between the rightmost local absorption maximum and the re-emission maximum.
However, there are exceptions to Kasha's rule, leading to the need for re-radiation matrices (Figure 2b), which characterize, for each incident wavelength (λ i ), a distribution of re-emission wavelengths In the context of path tracing, a normalized re-emission spectrum (column excluding λ i = λ o ) can be seen as a PDF of the wavelength shift of an irradiation λ i and a normalized absorption spectrum (row excluding λ i = λ o ) can be seen as a PDF of absorption of a wavelength λ i for a given outgoing radiation λ o .
To preserve energy conservation, the sum of each column of the re-radiation matrix never exceeds 1. The resulting value corresponds to the total albedo of the material at a given incident wavelength, that is, the sum of surface reflectance and fluorescence re-emission.
But, while re-radiation matrix is more flexible and suited for computer graphics applications, it has a significant memory footprint.

Gaussian mixture model
A GMM is a parametric probabilistic model representing a dataset as a linear superposition of Gaussian distributions. This linear superposition is represented by where: is the density of mixture of Gaussians, is a single Gaussian density, π k is the mixing coefficient (weight) of each Gaussian, K is the number of Gaussians. In this article, we model re-radiation matrices as mixtures of twodimensional multivariate Gaussian distributions: where: μ is the mean, is the covariance matrix.

Method Overview
Our technique consists of two related but distinct operations that form an efficient fluorescence rendering pipeline ( Figure 3): 1. Fitting, where the re-radiation matrices are fitted to GMMs (5). 2. Rendering, where the GMM representation is used to importance sample wavelength-shifting events (6).
The pipeline starts with a fitting preprocessing step. Here, the reradiation matrices representing fluorophores are filtered to remove invalid values and fitted to GMMs. During the rendering phase, we use these GMM representations for fluorescent wavelength-shifting events, where a conditional GMM is importance sampled. Depending on the path direction, we either sample based on absorption or re-emission wavelengths.
Using a GMM to represent a bispectral reflectance drastically reduces its memory footprint in the renderer, as instead of relying on a tabulated re-radiation matrix, we only need to store a few parameters of the mixture. This reduction is even more pronounced when the rendering technique relies on random sampling, as additional tabulated CDFs are necessary for efficient sampling of re-radiation matrices.

Fitting Phase
In the fitting phase, we apply parametric density estimation to fit a re-radiation matrix into a GMM. Later, we show that we can accurately re-construct the re-radiation matrix information from the GMM on the fly.

Dataset filtering
Some measurements in the dataset include a significant amount of noise, potentially lowering the accuracy of the fit. To mitigate this problem, we filtered the re-radiation matrices before the fitting phase.
Although not ideal, we zeroed values under the diagonal (λ i > λ o ) and negative values above it, since these values in our dataset correspond to measurement errors. We do not apply any further data filtering during the fitting phase. A further improvement would be possible with a proper characterization of the acquisition device's noise ratio.

Fitting algorithms
In our previous publication, we used the implementation of weighted EM provided by the Pomegranate library [Sch22]. In this revision, we instead use weighted EM in the better-behaving scikitlearn [PVG*11, sci21] library and compare its fitting accuracy with the weighted Bayesian fitting algorithm.

Weighted expectation-maximization
EM is a two-stage iterative algorithm suitable for the finding of mixture model parameters. The E-step (expectation) performs an estimation of the expected log-likelihood for a complete dataset. The M-step (maximization) maximizes the expected complete loglikelihood. The E-step and M-step are run iteratively until the expectation converges to the target. As a detailed introduction of EM is beyond the scope of our work, we review only a sketch of EM. For a detailed study on EM, please refer to Bishop [Bis06]. Vorba et al. [VKŠ*14] also showcase an adaptation of EM to the rendering context.
The EM algorithm works with a set of observation samples and computes a mixture matching its distribution. But, the re-radiation matrix instead provides weights of samples of pairs (λ i , λ o ). We could use a large set of samples distributed according to these weights, but then the accuracy of the representation would depend on the size of the sample set, whereas bigger sample sets would lead to higher memory consumption and a slower fitting process.
A more suitable approach is to use the weighted EM algorithm [GAPFH16, VKŠ*14] where all non-zero samples are used as observation outcomes and weighted according to their respective values.

Variational Bayesian inference
Variational inference follows the principle of the EM algorithm. It is a two-stage iterative algorithm that maximizes a posterior probability (MAP) instead of maximizing local likelihood (ML) as in the EM. Apart from updating the approximations in EM, the model evidence in variational inference also includes a pre-defined lowest prior distribution. By maximizing their posteriors using Bayesian inference in each iteration, the variational inference model gives more stable results compared to the EM, which can produce singularities and exhibit overfitting [Bis06]. Implementation of weighted variational Bayesian inference in our pipeline directly follows the derivation of weighted EM [HFW21].

Optimal number of Gaussians
Choosing an optimal number of Gaussians for the EM is a nontrivial problem, as with a higher number of Gaussians, we may face overfitting issues. Bayesian information criterion (BIC) can mitigate overfitting by the introduction of a penalty for the number of parameters of the model.
Variational Bayesian inference can take advantage of the predefined prior distribution, such as the Dirichlet process (used in the scikit-learn implementation as default) to select the optimal number of clusters based on the fitted data.
We evaluate our fitting on different numbers of Gaussians for EM as the selection may depend on the use case. Conversely, we let Bayesian inference use a lower number of Gaussians than the specified maximum if it opts to do so for improved accuracy of the fit.

Re-radiation re-construction
We can re-construct the re-radiation matrix from the GMM found previously. Gaussian mixture gives a probability density function, and so its integral is unity. Because the re-radiation integral does not hold this property, the resulting Gaussian mixture has to be scaled to re-construct the re-radiation later on: where: S is the scaling factor, (λ i , λ o ) is the re-radiation function, p(λ i , λ o ) is the probability density function of the GMM defined in Equation (1), are the wavelengths of a given absorption and re-emission event.
There are different strategies to retrieve the scaling factor S. We originally proposed two approaches:

Error minimisation
The first strategy consists of a minimisation process. S is set as a parameter to minimise the error between the measured re-radiation and the re-constructed mixture.
The choice of a norm function has little influence on the result when dealing with fluorescence -the data have a low dynamic range. By definition, this strategy reduces the average reconstruction error.

Integral ratio
The second strategy consists of defining S as the ratio between the integral of the original dataset and the GMM.
This method ensures that the re-constructed re-radiation has the same albedo as the original measured one.
Although this technique does not reduce the measurable reconstruction error, its application provides better rendering results because our eyes are more sensitive to brightness variation than to a slight chromaticity shift.

Strategy selection
The strategy for scaling factor computation depends on the target application. The integral ratio is better suited for perceptual applications, such as rendering with tone mapping. If the accuracy of the spectral re-construction of the re-radiation is needed, then error minimisation is a better choice.
In this publication, we apply integral ratio scaling as it is more closely related to the intended application of our method.

Rendering Phase
When a fluorescent event occurs during the rendering phase, we first need to sample an in-shifting or an out-shifting wavelength depending on the ray direction. Then, we need to evaluate the attenuation given a set of λ i and λ o wavelengths. The sampling phase relies on conditional probability: either λ i or λ o is known, and we have to sample a λ o or a λ i according to the direction of the ray, that is, evaluate an in-or out-shifting event.

Conditional GMM
To sample a random wavelength shift after a fluorescent event, we need to 'sample a slice' on the multivariate Gaussian mixture. With a fixed input wavelength λ i , we construct the conditional parameters for the resulting Gaussian mixture and generate the random variable as the output wavelength λ o .
Our model is made of: where μ k and k are the GMM parameters of the k th Gaussian.
We define T, the inversion of the kth covariance matrix as: The conditional probability that, given λ i , we shift to λ o in component k is: Then, the probability for the entire conditional GMM from a given λ i to λ o is: Finally, we generate the random variable based on the new conditional GMM model. First, we choose which Gaussian component we will sample using the conditional mixing coefficients π (λ i ) k , and then we generate the sample following the normal distribution of the chosen k th Gaussian: We only discuss sampling on the light path, but the same approach can be symmetrically applied on an eye path.

Importance sampling
When a ray hits a fluorophore, two events can occur: • the ray is reflected with the same wavelength as the incident ray -no fluorescent event, • the ray is absorbed and its energy re-radiated to a different wavelength -fluorescent event.
We first determine if a fluorescent event occurs based on the ratio between reflectance and fluorescence for a given incident wavelength. If a fluorescent event occurs, we use importance sampling to select a Gaussian distribution from the mixture and a specific wavelength from the distribution (Algorithm 1). Figure 4 shows an example of a conditional PDF of outgoing wavelengths of a specific incident wavelength λ i . Given an incident wavelength, we can generate this distribution from GMM and use it to importance sample the re-emitted wavelength λ o .
The probability of transition between absorption at wavelength λ i and re-emission at wavelength λ o is where: is the total reflected and re-emitted energy given λ i . Algorithm 1. Random wavelength shifting. Given a λ i , we either have no fluorescent event, that is, λ i = λ o or, we sample a λ o after a fluorescent event. In our renderer, we use HWSS [WND*14]. The HWSS technique multiplexes multiple wavelengths in a MC sample. It differs from a simple multiplexing approach by using a 'Hero' component, which is used to make all directional decisions. This technique drastically reduces colour noise in contexts with a wavelength dependence on the directional decisions, such as in a participating medium. In the case of fluorescence, we allow each of the wavelengths in the HWSS vector to be independently importance sampled. If done directly, we can find non-zero samples with a zero probability, which is incorrect for a MC integrator. Mojzík et al. [MFW18] introduced a balance term to avoid this artefact. We use the same term in our implementation: where: is the probability of jth HWSS sample, is the size of the HWSS vector.

Implementation details
To evaluate the re-radiation, we store the mean, covariance matrix and weight for each Gaussian of the mixture. We also need to store the scaling factor S. The determinant and inverse of the covariance matrix are pre-computed to speed up the evaluation of the mixture during rendering. The diagonal is stored in its original tabulated form without any alteration.
To efficiently perform conditional wavelength shifting, we also pre-compute two additional tabulated values to retrieve π (λ i ) k or π (λo) k : the sum of rows and the sum of columns of the Gaussian mixture. In our implementation, we use a 1-nm sampling rate.
Although additional pre-computed values increase the total memory footprint during rendering, it is still much smaller than the memory required when using tabulated CDFs, which are necessary for the use of re-radiation matrices.
We use rejection sampling to guarantee that absorption wavelengths do not exceed re-emission wavelengths. It proved sufficient in our case, as we observed only a small number of rejected samples. Alternatively, truncated distributions could be used.
The implementation of our method is available on Gitlab.

Results and Discussion
Our fluorophore representation consists of the original re-radiation matrix diagonal representing the reflectance spectrum and the scaling factor S with the GMM representing the fluorescent effect. We need to store seven values for each Gaussian of the mixture: • a mixing coefficient π (1 element), • a mean vector μ (2 elements), • a covariance matrix (2×2 elements).

Dataset
In our previous publication [HFW21], we evaluated our method on three measured re-radiation matrices distributed in the assets of the ART rendering toolkit -originally provided by Labsphere Inc. In this revision, we provide a more in-depth evaluation on a much larger bispectral dataset. We additionally use the dataset from RIT [GF00], which contains 144 measured re-radiation matrices of different materials.
All measurements in the resulting dataset come from the same device: Labsphere BFC-450. The materials were sampled at 10nm precision with an absorption wavelength ranging from 300 to 780 nm and a re-emission ranging from 380 to 780 nm, resulting in 49 × 41 re-radiation matrices.
One material measurement (IXCAXORA) was discarded from our evaluations. Its re-radiation matrix appears to be incorrect due to human error. We still include this sample in the supplemental.
In this section, we provide the results of the evaluation on the whole dataset but directly showcase only its small sub-set. Detailed results of our method on all samples are available in the supplemental.

Illuminants
To evaluate the accuracy of our representation, we compare the results of our fitted model with the original tabulated data. We base our accuracy evaluation on the CIE 2000 Delta E ( E * 00 ) between a rendering of the reference data and the corresponding fits observed under a large set of illuminants.
We use two sets of illuminants and distinguish the results of each group: • Monochromatic illuminants: We evaluate the bispectral reflectance of each fluorophore under monochromatic illumination from 300 to 780 nm to cover the full range of incident wavelengths in the dataset. The performance of our model is calculated as an average E * 00 of all monochromatic illuminants. • CIE standard illuminants: We also evaluate the colour reproduction accuracy of our model under the CIE standard illuminants, which cover often used illumination conditions in rendering practice.
Any measured differences between our model and pure reradiation matrices are caused by our compact representation of the fluorescence, as the reflectance spectrum is not modified.

Rendering and colour reproduction
While a MC rendering provides additional information about interreflection between different fluorescent materials and self-interreflection, the choice of the scene and its setup also limit the evaluation to particular configurations. And so, given the large number of material measurements, we instead predominantly use an analytical rendering method. Specifically, we evaluate the accuracy of the fitting on an analytically calculated interaction of light with a fluorescent material represented as a re-radiation matrix: where: C is the computed outgoing colour, λ i , λ o are the incident and outgoing wavelengths, is the colour matching function at λ o .
This equation is equivalent to a single light bounce configuration on a Lambertian surface where both incident and outgoing directions are aligned with the surface normal. This configuration allows us to assess the quality of the fits without any additional parameters interfering with our evaluations.
The re-radiation matrix corresponding to our fit can be constructed by evaluating the GMM at regular intervals for λ i , λ o . This approach allows efficient generation of noise-free results for each material in the dataset for a large variety of illumination types.

Rendering accuracy
Due to the fluorescent nature of our dataset, illumination conditions have a significant impact on error, making an objective evaluation difficult. Therefore, we evaluate our representation of fluorophores under a wide variety of illuminants visualized as • Colour band: To visualize outgoing colour after interaction of a monochromatic illumination with a fluorophore, we use a 'colour band'. It consists of three strips representing from top to bottom: outgoing colour based on a re-radiation matrix, outgoing colour based on our fitted model, and visualization of E * 00 between these two. The colour bands cover the entire excitation spectrum of the fluorophore (horizontal axis) and allow wavelengthspecific comparisons of the measured bispectral reflectance with our fit. • Colour patches: To compare the quality of our fit with the original measured data, we use colour patches to visualize how well a fluorescent material is represented under a specific CIE standard illuminant. The fitted patch is embedded inside of a reference patch. The colour values are divided by the Y value of the used illuminant. Figure 5 shows the five worst fluorophore fits based on their average E * 00 for monochromatic illuminants at four Gaussians using EM. While the colour difference is visible for some monochromatic illuminants, the difference is hardly noticeable for most CIE standard illuminants, especially with a higher number of Gaussians. While four Gaussians would be suitable for most rendering purposes, eight Gaussians already provide high accuracy within a small memory footprint. This also applies to global illumination, as showcased in Figure 1. This is supported by a comparison of the quality of the GMM fit with respect to the number of Gaussians used ( Figure 6). Even though the reconstruction of the original re-radiation at four and two Gaussians is very rough, it does not translate to big observable differences as our visual system is much more sensitive to luminance discrepancies than to chromaticity variations.

Average accuracy on the whole dataset
To compare the overall accuracy of the two fitting methods with a varying number of Gaussians, we provide a E * 00 averaged across the entire dataset. Figures 7b and 8b show a consistent decrease of the E * 00 with an increasing number of Gaussians. Bayesian inference performs better with a lower number of Gaussians, but is outperformed by EM from eight Gaussians. This is caused by the Bayesian inference implementation opting to use a lower number of Gaussians for fitting than the allowed maximum. Figure 7a shows the average E * 00 over the entire dataset for monochromatic illuminants ranging from 300 to 780 nm. The significant error observed in the lower parts of the spectrum can be attributed to limitations of the E * 00 involving colours with low luminance. In our case, this part of the spectrum usually represents non-visible light or a negligible amount of re-radiation, as absorption bands of fluorophores are often situated there. Therefore, the practical influence of this error on real-world rendering applications is limited. The monochromatic illumination error is visualized for a few selected fluorophores in Figure 5 and the whole dataset in the supplemental.

Model accuracy
In Figure 8a, we provide an evaluation of a more practical application of our approach-rendering with CIE standard illuminants. We show an average of the whole dataset for the ten worst performing illuminants while fitting four Gaussians with EM. On average, we achieve a E * 00 lower than 1.1 for two, 0.6 for four and 0.3 for eight Gaussians.

Memory footprint
While retaining high accuracy, our method is memory efficient, requiring only a fraction of the memory necessary for use of the raw re-radiation matrices.
In the most favourable case, we only have to store the upper triangular part of the re-radiation matrix due to the nature of fluorescent re-radiation being of lower energy than excitation caused by the incident illumination. For our dataset measured at 10 nm precision, it is approximately 1150 values. If maximal precision was required, each matrix could be interpolated to 1 nm precision, resulting in a hundredfold increase in memory footprint. Additionally, to facilitate efficient importance sampling of wavelength-shifting events, it is desirable to compute tabulated CDFs, effectively tripling the total required memory when using re-radiation matrices for rendering.  With our method, we only need a small set of GMM parameters to both evaluate the re-radiation on the fly and to importance sample wavelength shifting events.
In practice, for a representation using N Gaussians, we need to store 1 + 7 · N values. For importance sampling, we need two additional 1D arrays to efficiently evaluate the conditional GMM. This is much less than the two tabulated 2D CDFs required for a fluorescent material when using a re-radiation matrix. Additionally, if desired, our method can be used solely for importance sampling instead of tabulated CDFs.
The memory efficiency of our method becomes especially pronounced when working with larger numbers of re-radiation  matrices. It makes fluorescent textures where each pixel may have to be represented by a re-radiation matrix usable in rendering practice. This, in turn, facilitates the use of related techniques such as enlargement of the colour gamut by uplifting of RGB textures [JWH*19].

Conclusion and Future Work
We presented a method for efficient storage of re-radiation matrices via GMM. We showed that this approach drastically lowers required memory for the rendering of fluorescence while retaining colour reproduction accuracy matching that of the original data. Further, the resulting GMM representation can be directly used in MC-based rendering pipelines, as it supports efficient importance sampling.
We evaluated weighted EM and Bayesian inference fitting approaches on an extensive dataset of measured re-radiation matrices under all CIE standard illuminants and a wide range of monochromatic illuminants. Both fitting methods proved to provide a comparable and highly accurate representation of fluorescent reradiation with as little as eight Gaussians for most practical rendering use cases.
Different models and fitting approaches could be explored, but considering the rendering accuracy and small memory footprint of the presented solution, the potential gains seem marginal.
We think that the main direction of future work is in the development of a parametric fluorescence model. Such model would not require tabulated CDFs -further lowering total memory requirements for efficient rendering of fluorophores. More importantly, it would allow a smooth integration of editable fluorescent materials into artistic workflows.