1. Introduction

Monte Carlo simulations of light transport in natural materials are usually performed under the assumption that such materials may be represented using layers and that the directional changes of travelling rays can be modeled using precomputed phase functions. In this paper, we provide the mathematical framework for an alternative approach in which these directional changes are computed on the fly, without explicitly storing the material constituents. This approach strikes a balance between the use of precomputed phase functions [1, 2], which do not relate to any particular material, and a conventional ray tracing approach, which, among other difficulties, requires much storage space if one wants to geometrically represent the material constituents [3]. Although in the context of this paper we focus on its application to the modeling of sand optical properties, this approach can potentially be applied to other materials.

In this paper, we present a spectral light transport model for sand, hereafter referred to as SPLITS. We evaluate the model using virtual spectrophotometric [4] and virtual goniophotometric [5] methods, and comparing its results to measured data. SPLITS represents, to the best of our knowledge, the first attempt to simulate the spectral reflection properties of sand using its physical and mineralogical properties, henceforth referred to as characterization data, as input.

There are several reasons for modeling sand in terms of its characterization data rather than in terms of arbitrary parameters not directly relating to sand. Hyperspectral remote sensing with satellite or aircraft based equipment is used to investigate properties of land surfaces without having to physically survey the area [6, 7]. For other planetary bodies, such as Mars, a field survey may not be possible. In this case, spectral signatures may be the only clue available for determining surface characteristics [8, 9]. For these scientific applications, this relationship between model parameters and the physical properties is necessary so that sand characterization data may be derived from predictive simulations provided by the model. Finally, for the model to be evaluated objectively, the parameters must relate to the material that one is trying to simulate.

The remainder of this paper is organized as follows. Background information about sand is given in the next section. In Section 3, previous work relating to the modeling of sand optical properties is outlined. The model is described in Section 4. The evaluation approach and the results are discussed in Section 5. We conclude and outline directions for future work in Section 6.

2. Properties of Sand

In this section, the relevant background information on sand is provided. We begin by defining precisely what sand is. The pertinent factors affecting light transport with sand are then described. The following, however, is not a complete treatise on sand. Interested readers are referred to works by Pettijohn et al. [10], Brady [11], or Gerrard [12].

Sand is a particular type of soil. Soil is composed of particles of weathered rock and sometimes organic matter immersed in a medium of air and water (the pore space) [12].

Soils are classified according to the size distribution of the mineral particles [11]. This is accomplished first by assigning individual particles to classes, called soil separates, according to their size. The relative masses of each of the soil separates are then compared to determine the texture of a soil sample.

Various agencies have differing definitions for soil separates and textural classes [11]. In this work, we use the system developed by the United States Department of Agriculture (USDA) [13]. The USDA defines three soil separates, called sand, silt, and clay, as delineated in Table 1. Particles larger than 2mm are classified as rock fragments and are not considered to be part of the soil. A sand textured soil contains at least 85% sand-sized particles.

That the term sand is used to describe both a soil separate and a soil texture may be a source of confusion. In the remainder of this paper, the term sand is used to refer to a soil texture unless otherwise stated (for instance, by referring to sand-sized particles).

Table 1:. Soil separates (particle size classes) defined by the United States Department of Agriculture [13].

View Table | View all tables in this article

2.1. Factors Affecting Light Transport

The reflectance of sand generally increases with wavelength [14]. Additionally, sand reflection is non-Lambertian [15, 16]. In particular, sand exhibits retro-reflection [16], which is reflection in the direction toward the source of illumination. Sand also reflects light in the forward direction, peaking at an angle to the normal greater than that of the direction of specular reflection [16].

The most important factors contributing to the reflectance of soils in the visible range are its mineral composition (and iron oxides in particular), moisture, and particle size and shape [17].

2.1.1. Mineral Composition

As sand is composed primarily of weathered rock [10], the optical properties of that rock may influence light transport within sand. The parent material is the rock that is the source of the mineral part of the sand. This is typically a silicate mineral [11] such as quartz, gypsum or calcite, with quartz being the most common [11, 18]. While these minerals are colourless in pure form, trace amounts of contaminants may substantially affect their colour [19].

Iron oxide gives sand its distinctive hues. Indeed, the two most significant minerals that determine soil colour, hematite and goethite [20], are iron oxides. Goethite, also known as yellow ochre [20] or limonite [21], is one of the most common minerals found in soils [20]. It colours soils yellow to brown [22]. Hematite, or red ochre [20], imparts a red colour to soils and may mask the colour of goethite except when in small quantities [22]. Hematite is usually found in tropical regions and is rare in temperate or cool climates [23]. Iron oxides may be present as contaminants in the parent material [21]. They are also found, typically within a kaolinite or illite matrix, as coatings, approximately 1–5µm thick, that form on the grains during aeolian (i.e., by wind) transport [24].

Additionally, magnetite and ilmenite are often found in beach and river sands [25]. These minerals are spectrally similar. They are opaque and are black in colour [25].

2.1.2. Water

The presence of water darkens sand. The principal reason for this is that the reduced contrast between the refractive index of the pore space and that of the mineral particles (typically quartz) induces stronger forward scattering at particle interfaces [26, 27].

2.1.3. Grain Size and Shape

Grain size affects reflectance by influencing the number of scattering interfaces per unit distance through the medium [28]. Smaller particles, and thus a higher density of scattering interfaces, result in higher reflectances [17, 18].

Grain shape may also affect scattering properties [17]. The shape of the grain is described by two quantities: sphericity and roundness. Sphericity refers to the general shape of a particle by expressing its similarity to that of a sphere [29]. In this paper, we use the sphericity measure proposed by Riley [30], which is a projection sphericity measure, meaning its definition is based on the projection of the particle onto a plane. The Riley sphericity, Ψ, of a particle is given by

where D_i is the diameter of the largest inscribed circle and D_c is the diameter of the smallest circumscribed circle, as shown in Fig. 1.

 

Fig. 1. Projection of a sand grain onto a plane. The circles indicate measurements used in the calculation of the Riley sphericity [30].

Download Full Size | PDF

In contrast, roundness can loosely be described as a measure of detail in the features on the grain surface [29]. A higher roundness value indicates a smoother surface. Roundness is determined by comparing the particle to images on a roundness chart [31].

2.1.4. Additional Factors

There are several other factors that may effect the reflectance of soils, such as organic matter content, surface roughness due to aggregation of soil particles, or human activities like tillage or transportation. These factors, however, either show themselves primarily outside the visible region of the spectrum, are less applicable to sands than to finer soils, or affect reflectance indirectly.

3. Related Work

In this section, we provide an overview of general purpose models which simulate optical characteristics relating to the appearance of sand. We then outline methods available for simulating some of the factors affecting light transport in sand described in Section 2.1. The reader interested in more information about general particulate scattering models based on radiative transfer theory is referred to the recent works by Zhang and Voss [32], Xie et al. [33], and Mishchenko et al. [34, 35].

3.1. General Models

Hapke [36] proposed an analytic radiative transfer model for particulate surfaces, which was then used to simulate reflectances from planetary surfaces [37]. Emslie and Aronson [38] modeled particulate surfaces, treating large particles separately from those smaller than the wavelength of light. Large particles were modeled using spheres of the same volume as actual particles and absorption was simulated at small scale asperities on the surface of the particle. However, their model was only applicable to the far infrared [39]. Egan and Hilgeman [39] subsequently rectified this by simulating scattering, in addition to absorption, at these asperities. Oren and Nayar [40] generalized the widely used Lambertian model [41] by representing a surface using a collection of symmetrical V-shaped cavities. Oren and Nayar compared the output from their model to a sand sample. Their model, however, only simulates reflectance in the spatial domain. Mishchenko et al. [42] proposed a radiative transfer technique for computing the bidirectional reflectance of a semi-infinite random media composed of nonabsorbing or weakly absorbing particles. Stankevich and Shkuratov [43] simulate light scattering in optically inhomogeneous particulate media using a geometric optics approach which combines ray tracing techniques and Monte Carlo methods. Recently, Peltoniemi [44] used a similar framework to simulate light scattering in densely packed particulate media and account for polarization. It is worth noting that the preceding models are intended to simulate a wide variety of materials. As such, their parameters do not specifically relate to sand.

3.2. Simulation of Specific Effects

3.2.1. Moisture

Various methods exist to achieve the darkening effect caused by water. Jensen et al. [45] use an extension of the Henyey-Greenstein phase function [46], and adjust the degree of forward scattering to achieve varying levels of wetness. We remark that the Henyey-Greenstein phase function is an empirical function with a single parameter which bears no relation to the sand characterization data [47]. It is also worth noting that it has been demonstrated by Mishchenko et al. [42] that its application in the simulation of soil scattering properties results in large errors, which can exceed a factor of twenty at backscattering geometries and a factor of three near forward scattering geometries. Lobell and Asner [48] use an exponential function to empirically model reflectance in terms of moisture content. Neema et al. [49] analytically model the sand medium using layers of reflective spheres coated by a thin film of water, the thickness of which is a function of moisture content.

3.2.2. Iron Oxides

Barron and Montealegre [50] apply Kubelka-Munk theory [51] to model the influence of iron oxides on soil colour. Torrent et al. [22] relate hematite and goethite concentrations to a function of the Munsell colour [52] of the soil, which Torrent et al. call the redness rating. Okin and Painter [53] use spherical, hematite coated quartz particles to study the effect of grain size on the reflectance of desert sands.

4. The SPLITS Model

The SPLITS model is a comprehensive light transport model for sand. That is, both the spectral (colour) and spatial (bidirectional reflectance distribution function, or BRDF [54]) responses are simulated. These two components together constitute the measurement of appearance [55]. Furthermore, the SPLITS model is dependent on the physical and mineralogical characterization data for sand. The appearance of sand depends on many such parameters, and, to the best of our knowledge, there is no single source describing all of the requied characterization data for several sand samples. Therefore, data was gathered from several sources. With any research effort, this data gathering process represents a cruicial step. As pointed out by several researchers, “good science requires both theory and data — one is of little use without the other” [56].

The remainder of this section is organized as follows. We begin with a general overview of the model in Section 4.1. The construction of the model is described in Section 4.2. In Section 4.3, the light transport simulation is described. We conclude with an outline of possible extensions supported by the model in Section 4.4.

 

Fig. 2. Modeled sand medium.

Download Full Size | PDF

4.1. Concept

The SPLITS model consists of sand particles that are randomly distributed within the half space below a horizontal plane representing the sand surface. The pore space, defined in Section 2, consists of a mixture of air and water. This is depicted in Fig. 2. The particles are composed of the most important mineral constituents of sands: quartz, hematite, goethite and magnetite. In addition, quartz particles may be contaminated by hematite or goethite, or coated by hematite or goethite in a kaolinite matrix (Section 2.1.1).

Light propagation in the model is described in terms of ray (geometric) optics. Its implementation is based on an algorithmic process combining ray tracing techniques and standard Monte Carlo methods, which are applied to generate a scattered ray given the incident ray and the properties (e.g., wavelength dependent refractive indices) of the sand medium (Fig. 2). The wavelength of light, a wave (physical) optics parameter, is included in its formulation by associating a wavelength value with each ray. It is assumed that each ray carries the same amount of radiant power, and the energies associated with different wavelengths are decoupled (i.e., phenomena such as fluorescence are not addressed). Although the geometric approach adopted in the model design does not favour the direct computation of wave optics phenomena, such as interference¹, it is more intuitive and allows an adequate description of the large scale light behaviours, such as reflection and refraction, in environments characterized by incoherent radiation fields, which represent the main targets of our simulations. It is also assumed that the relevant distances are much larger than the wavelength of the light. This assumption holds over the visible region of the spectrum for sand particles, as they are larger than 0.05mmin diameter (Section 2).

Ultimately, an exact simulation of a particulate material would involve computing the solution of Maxwell’s equations [58]. However, this is not feasible given the current state of computer technology and the complexity of the media involved. Any model must therefore be approximate: either in that the medium being simulated is an ideal one, or in that the light transport simulation is approximate. For instance, it could be argued that an approach based on Mie theory may yield better results. However, simulating Mie theory may add undue complexity and, due to the large size and irregular shape of sand particles, may not result in an improvement over a ray optics approximation [59].

Table 2. A summary of the data required in the simulations, along with typical ranges and values used.

View Table | View all tables in this article

4.2. Construction

We begin by describing the representation of the medium in which the particles are immersed. Next, details are provided for the geometry and composition of individual particles. We then describe how those particles are distributed by size, shape, and composition. Table 2 lists the data used by the model along with typical ranges and the values used.

4.2.1. Extended Boundaries

As previously stated, the sand medium is bounded above by a horizontal plane representing the surface. The medium is additionally bounded below by another horizontal plane to prevent runaway paths from occurring during the light transport simulation. These two planes are called the extended boundaries (Fig. 2). The distance, D, between them is set high enough so that light penetration to that depth is negligible [49].

4.2.2. Pore Space

The pore space, as previously stated, consists of the air and water between the sand particles. The porosity, P, defines the volume of pore space as a fraction of sand volume. The volume of water in the sand is expressed as the fraction, S, of the porosity, and is called the degree of saturation. The amount of water in the sand may also be expressed as a volumetric water content, w, which is the fraction of the total volume occupied by water [48]. This may be related to the degree of saturation by S=w/P.

 

Fig. 3. An example of the light propagation through a modeled sand particle. Left: The particle consists of a core and an optional coating. The small box around the ray path at the coating interface corresponds to the layer interface diagram to the right. Right: A closer look at the modeled interface between the core, coating, and/or the surrounding medium.

Download Full Size | PDF

4.2.3. Particle Geometry

Individual particles consist of a core and an optional coating, as depicted in Fig. 3 (Left).

The core is modeled as a prolate spheroid with the semiminor axis of length a and the semi-major axis of length c, with 0<a≤c, as depicted in Fig. 4. The particle diameter is thus given by s=2c. When placed on a horizontal plane, a prolate spheroid will naturally lay so that the major axis is parallel to the plane. The projection of the prolate spheroid onto a horizontal plane is therefore an ellipse with semiminor and semimajor axes of a and c respectively. Thus, from Eq. (1), a and c are related to the Riley sphericity, Ψ, by the expression

One might suggest that because sand particles are irregularly shaped, one should use more complex particle shapes. To justify the added model complexity, it would be necessary to have the data describing the particle shapes at that level of detail. To the best of our knowledge, such data is not available in the literature. Since roundness and Riley sphericity data is readily available [65], we decided to use prolate spheroids. However, SPLITS is not limited to the use of prolate spheroids. The use of arbitrary particle distributions with SPLITS is described by Kimmel [66]. Optionally, the core may be coated by a layer of uniform thickness, h, that is proportional to the particle size, with h′=h/2c. The coating thickness is assumed to be small relative to the size of the particle, so that light travelling within the coating does not stray far from the point of entry. Hence, the coating may be approximated locally as a flat slab.

The interfaces between the core, coating, and surrounding medium are modeled using randomly oriented microfacets of equal area to simulate a rough surface, as shown in Fig. 3 (Right). The orientations of the facets are distributed such that the dot product between the microfacet normal, n′, and the interface normal, n, is given by

 

Fig. 4. Prolate spheroid.

Download Full Size | PDF

where X is normally distributed with zero mean and standard deviation

$\sigma =\frac{1-R}{2}.$

This standard deviation was chosen so that R<n′·n≤1 for 95% of the facets [67]. Additionally, the microfacet normals are constrained so that n′·n>0. The particle roundness, R, is therefore used to control roughness. When R=1, the interface reduces to a smooth surface, as one would expect based on the concept of roundness [29].

4.2.4. Particle Composition

The minerals that are used in the model are quartz, hematite, goethite, and magnetite. Additionally, kaolinite is used as the coating matrix (Section 2.1.1). These minerals may occur in pure form or as a mixture.

To represent a mixed material, we use an equation from Maxwell Garnett theory [68] that relates the dielectric constant, the square of the complex refractive index, of a mixture to those of its constituents. This theory originally addressed the optical properties of media containing minute metal spheres and it was developed to predict the colours of metal glasses and metallic films [69]. Recently, it was applied to assess the effect of grain size on the reflectance of sandy materials [53]. According to Maxwell Garnett theory, the dielectric constant, ε_avg, of a mixture is given by

where ε_m is the dielectric constant of the matrix, ε is the dielectric constant of the inclusions, and ν_i is the volume fraction of the inclusions [68].

4.2.5. Particle Types

The sand particles are divided among the following eight types: pure (quartz (pq), hematite (ph), goethite (pg), and magnetite (pm)), mixed (hematite with quartz (mh) and goethite with quartz (mg)), and coated (hematite coated quartz (ch) and goethite coated quartz (cg)). The hematite and goethite in the coatings are present in the form of inclusions in a kaolinite matrix (Section 2.1.1).

Ultimately, we must determine the fraction of volume, ν_j(1-P), occupied by particles of type j, as j varies over each of the eight particle types (pq, ph, pg, pm, mh, mg, ch, cg). For the mixed particles, we must also determine the volume concentration, ν_j,ħ, of each of the constituent minerals ħ within the particle, required to apply Eq. (4). For the coated particles, we must determine the volume concentration of the inclusions within the coating, which may be computed from the overall volume concentrations, ν_j,ħ, of each of the mineral constituents of the particle. We begin by determining the mass fraction, µj, of each of the various types of particles (Table 3), and the mass concentration, ϑ_j,ħ, of each mineral, ħ, within particles of each type, j (Table 4).

The overall mass concentrations of hematite, goethite, and magnetite in the simulated sand is controlled, respectively, by ϑ_h, ϑ_g, and ϑ_m. In addition, we define some additional quantities. We define the iron oxide concentration as

Because of the particular importance of hematite and goethite (Section 2.1.1), we also define the concentration of these two minerals,

and a ratio describing the relative proportions of hematite and goethite,

The concentrations of hematite and goethite may be given equivalently by the above two quantities. For simplicity, the concentration of hematite (respectively goethite) in the mixed and coated particles containing hematite (respectively goethite) is fixed at ϑ_hg. The remainder of the mineral matter is quartz (within particle cores) and kaolinite (within coatings).

Three parameters, µ′_p, µ′_m, and µ′_c, partition the particles by mass into the pure, mixed, and coated particles respectively, with µ′_p+µ′_m+µ′_c=1. These parameters are further constrained by the concentrations of the various mineral constituents, since, for example, a particle consisting of a quartz core coated by a mixture of hematite and kaolinite has an upper bound on hematite concentration within that particle.

From the above parameters, we may now derive the mass fraction, µ_j, of each of the various categories, j, of particles. These are indicated in Table 3, with two unknowns, β_q and β_hg. Since we are given the total mass concentration of hematite, ϑ_h, and the mass concentration of hematite in each of the particle types (provided in Table 4), we have the relationship

r_hgβ_hg+ϑ_hgr_hgµ′_m+ϑ_hgr_hgµ′_c=ϑ_h.

From Equations (6) and (7), we see that ϑ_hgr_hg=ϑ_h. Substituting and solving for β_hg yields ${\beta }_{hg}=\frac{{\vartheta }_h(1-{\mu }_m^{\text{'}}-{\mu }_c^{\text{'}})}{r_{hg}}.$

Noting that µ′_p=1-µ′_m-µ′_c, we get

The reader may wish to verify that we arrive at the same result if we relate the concentration of goethite in each particle types to the total goethite concentration. Now we may solve for β_q. Since top row in Table 3 must sum to µ′_p, we have

Table 3. Mass fractions, µ_j, of each of the different classes, j, of particles. The quantities β_hg and β_q are unknowns, given by Eqs. (8) and (9) respectively.

View Table | View all tables in this article

Table 4. Mass concentration, ϑ_j,ℓ, of the constituent minerals, ℓ, in each of the different classes j of particles. The quantities ϑ′_k and ϑ″_k are unknowns given by Eqs. (10) and (11) respectively.

View Table | View all tables in this article

Next, we must determine the two unknowns, ϑ′_k and ϑ′_k, in Table 4. These quantities correspond to mass concentrations of kaolinite in the hematite coated quartz and the goethite coated quartz particles, respectively. Recall that the hematite coated quartz particles consist of a pure quartz core coated in a mixture of hematite and kaolinite. Hence, to determine ϑ′_k, corresponding to the mass concentration of kaolinite within the particle (Table 4), we need to know the volume of the coating as a fraction of the total volume of the particle, ν_coat. Solving for ϑ′_k [66] then yields

Similarly,

The volume of the coating may be estimated as A_Sh, where A_S is the surface area of the particle, since we are assuming that the h≪s. Therefore,

${\nu }_{\mathrm{coat}}\approx \frac{A_S(s,\Psi )h}{V(s,\Psi )+A_S(s,\Psi )h},$

where V is the volume of the particle. Dividing the numerator and denominator by V yields

where AV (s,Ψ)=A_S(s,Ψ)/V(s,Ψ) is the surface area to volume ratio of the particle. Noting that

$A_V(s,\Psi )=\frac{A_S(s,\Psi )}{V(s,\Psi )}=\frac{s^2{A}_S(1,\Psi )}{s^{3}V(1,\Psi )}=s^{-1}{A}_V(1,\Psi ),$

and that, by definition, h=h′s the particle size,s, may be eliminated from Eq. (12), yielding

The surface area to volume ratio of a prolate spheroid with s=2c=1 is given by [66]

For simplicity, we use the mean sphericity, $\overline{\Psi }$ , in the above (i.e., AV (1,$\overline{\Psi }$ )) rather than having ν_coat, and thus ϑ′_k and ϑ′_k, vary with sphericity.

Now that we know, for each particle type j, the mass concentration ϑ_j,ħ of each mineral ħ (Table 4), we may compute the density of a particle of type j,

Similarly, we may compute the particle density, γ (i.e., the mean density over all particles), knowing the mass fractions µ_j of each particle type j (Table 3). This is given by

Knowing the mass fractions and densities of the components of a mixture allows us to compute the volume fractions. Thus, we may now compute the volume concentrations of all the minerals within a particle, as well as the volume fractions of each particle type. The former is given by

and the latter by

For hematite coated particles, the volume fraction of hematite within the coating, required to compute the complex refractive index of the mixture of hematite and kaolinite, is ν_ch,_h/(ν_ch,_k+ν_ch,_h). The volume fraction of goethite within the coating of the goethite coated particles is determined similarly.

4.2.6. Shape and Size Distribution

The size and shape of the particles are randomly distributed and are independent of one another. That is, the conditional size distribution for any two shapes is the same.

The sphericity is normally distributed, with the mean and the standard deviation derived from data provided by Vepraskas and Cassel [65], and presented in Table 2. The sphericity is also constrained to fall within a range derived from the same data.

The particle size is distributed according to a piecewise (one piece for each soil separate) log-normal distribution as suggested by Shirazi et al. [70]. That is, logs is normally distributed. This distribution is characterized by two parameters: the geometric mean particle diameter, d _g, and its standard deviation σ_g, which are functions of soil texture (as defined in Section 2). It is important to note that, since the distribution is specified in a piecewise manner, there will be a different d_g and σ_g for each soil separate. Defining

where a_g=logd_g and b_g=logσ_g, the mass fraction of the particles with sizes ranging from s1 to s2 is then given by

Assuming that particle density does not vary significantly with size, 𝕱_m(s1, s2) may also be interpreted as a volume fraction. This approximation is justified by the fact that the silt and sand-sized particles are dominated by quartz [11]. Shirazi et al. [70] provide a table for determining d_g and σ_g from texture.

4.3. Light Transport Simulation

When a ray is incident at the outer extended boundary of the sand medium, it is either reflected or refracted, with the probability of reflection given by the Fresnel equations [58], ² according to the following algorithm. After computing the Fresnel coefficient at the boundary, a random number uniformly distributed on (0,1) is computed. If the Fresnel coefficient is greater than the random number, then a reflected ray is generated applying the law of reflection. Otherwise, a refracted ray is generated according to Snell’s Law.

Once a ray penetrates the outer boundary, it has entered the medium. The ray is then repeatedly scattered by the sand particles until it escapes the medium or is absorbed. Rather than explicitly represent and store each particle, however, they are generated stochastically as needed and subsequently discarded. The light transport simulation process is depicted in Fig. 5.

4.3.1. Extended Boundaries

A ray reaching the surface extended boundary from the inside, as from the outside, is reflected with a probability given by the Fresnel equations [58] as described above. If internally reflected, the ray continues traversing the medium. Otherwise, the ray escapes the medium and the scattered ray is returned. When a ray reaches the lower extended boundary, it is absorbed.

4.3.2. Pore Space

Each time a ray reaches an interface to the pore space, the medium representing the pore space is selected randomly. The degree of saturation, S, defined in Section 4.2.2, is used to partition the pore space into air and water. Water is selected with probability S to represent the pore space. Air is selected with probability 1-S. This selection is made when an incident ray approaches, as well as when the ray reaches the outer interface of a particle from the inside.

4.3.3. Generating a Sand Particle

As previously stated (Section 4.1), we do not explicitly store each sand particle. Rather, they are generated on the fly. The light interaction with this particle is then simulated, and the particle is then discarded. In a conventional ray tracing approach, particles would be stored explicitly. Random variables, such as the distance to the next particle, the shape, size, and composition of the particle that is intercepted, and the point on the surface of the particle that is intercepted, arise implicitly as a consequence of the ray tracing simulation. Conversely, in the SPLITS model, we compute the probability distribution functions for these random variables and sample them explicitly. In addition, we ensure that the particle lies completely between the extended boundaries and that the particle does not intersect with the previous ray.

 

Fig. 5. Flowchart outlining the light transport algorithm used by the SPLITS model. The dashed line indicates the main loop. The term inner boundary refers to the inside of either extended boundary, whereas the outer boundary is the outside of the surface extended boundary.

Download Full Size | PDF

The physical distance to the next particle is determined by randomly choosing a distance for each category j of particles and then selecting the category for which that distance is minimal. The distance, d _j, to the next particle of type j is an exponential random variable with a mean of 1/K_j, where K_j is the geometric attenuation coefficient for particles of type j. This distance is given by

where ξ is uniformly distributed on (0,1) [67].

For prolate spheroids with normally distributed sphericity and log-normally distributed size, K_j reduces [66] to

where

and

The particle size distribution, f_g(s), is given by Eq. (19), AV (1,Ψ) is given by Eq. (14), and Φ′(x̄,σ²)(x) is the probability density function for the normal distribution with mean x̄ and variance σ² [67].

The particle size, s, is then chosen according to the probability density function

and the sphericity, Ψ, is chosen according to

where C ₁ and C ₂ are the constants that ensure that f_s and fΨ (respectively) integrate to one over their domain [66]. Note that Eq. (25) does not match the distribution provided by Shirazi et al. [70] shown in the integrand of Eq. (20). This follows from the distinction between the particle size distribution by mass, provided by Shirazi et al. [70] (or by volume as we are interpreting it), and the distribution of particle sizes struck by rays travelling randomly through the medium [66].

Now that we have the particle size, s, the coating thickness is then given by h=h′s, where h′ is the relative coating thickness specified in Table 2. Additionally, the particle roundness, R is given by a normal random variable with mean R̄ and standard deviation σ_R given in Table 2.

The point on the surface of the particle which is struck is chosen uniformly on the side of the particle surface facing the ray origin. To select a point uniformly on a prolate spheroid (Fig. 4), we set z=cF ^-1(2ξ₁-1), where ξ₁ is a uniform random number on (0,1), F(u) is the fraction of the surface area of the spheroid above the plane z=u, given by

and $e=\sqrt{1-{\Psi }^4}$ is the eccentricity of spheroid [71]. Since F ^-1 is difficult to compute analytically, F may be inverted using numerical techniques. The other two coordinates are then given by x=ar cosθ and y=ar sinθ, with θ=2πξ ₂ (ξ₂ being another canonical random variable) and $r=\sqrt{1-z^2⁄c^2}$. To restrict the generated point to the side of the particle facing the ray origin, we transform the ray direction v into the particle’s local coordinate space to get v′, and compute the normal n at the chosen point. We then use the point (x,y,z) if n·v′<0 and (-x,-y,-z) otherwise.

Before we simulate light interaction with this particle, we first check the validity of the particle. If the point, q, that is d units along the ray lies outside the extended boundaries, the ray will instead interact with the boundary as described in Section 4.3.1. If the randomly generated particle intersects with either boundary, the particle is rejected. To account for the opposition effect [72] (i.e., the increase in reflectance toward the source of illumination due to shadow hiding), the particle is also rejected if it intersects with the last leg of the path. If the particle is rejected, the above process of generating a distance and particle intersection point is repeated.

4.3.4. Light Propagation Within a Sand Particle

Light propagation within a sand particle is simulated using a random walk process as exemplified in Fig. 3. A ray incident upon the particle first interacts with the outer interface: either between the pore space and the coating or between the pore space and the core if no coating is present. The ray may be reflected or transmitted at each interface. When the ray is transmitted into the core (or internally reflected back into the core) a ray-particle intersection is computed to determine the next point at which the ray may exit the core. Absorption is simulated within the coating and within the core. This process is repeated until the ray escapes the particle or is absorbed.

At each interface, the microfacet normal is selected according to Eq. (3). As mentioned in Section 4.2.3, the microfacets are of equal area, A. The projected area with respect to the incident ray, v, is therefore A|n′·v|. Hence, the probability that an incident ray v strikes a microfacet with normal n′ should be scaled by |n′·v|. This is accomplished using the rejection method [67]. The ray is then reflected in n′ with a probability once again determined using the Fresnel equations³ [58], considering the media on either side of the interface. Otherwise, the ray is refracted according to Snell’s Law for absorbing media [73]. If the ray was to be reflected and n·r<0 (where r is the reflected vector), or if the ray was to be refracted and n·t>0, then multiple scattering is approximated using a cosine lobe centered around n (respectively -n) [41].

Given the extinction index, k, of the coating or of the core, absorption within the coating or within the core is then simulated according to Lambert’s Law [54]. The light is transmitted a distance d with probability T=exp(-αd), where α is the absorption coefficient, given by

$\alpha =\frac{4\pi \Im \left(\eta \right)}{k}$

[58]. For the coating, d=h/|n·t|, where t is the ray transmitted through the coating. For the core, d is determined from the ray-particle intersection.

4.4. Extensibility

The framework for SPLITS supports several possible extensions. The model may be extended to include other minerals, given their densities and spectral complex refractive indices. The model can also be extended to include other particle shape and size distributions [66]. Additionally, multiple coatings may be applied to the particles.

While polarization was not simulated, it can be incorporated into the SPLITS model by tracing the Stokes vector instead of the intensity. It is worth mentioning that during the model development an effect of polarization was examined, namely birefringence, since it may affect the accuracy/cost ratio of the simulations. In the case of quartz, for example, the birefringence is low [74], and therefore only the ordinary ray refractive index was used. A full description of the data and optics equations used to implement the SPLITS model is provided by Kimmel [66].

5. Evaluation

To evaluate the SPLITS model objectively, it is necessary to perform comparisons with measured data directly. Thus, comparisons were made between measured data [75] and results obtained from the model using virtual goniophotometric [5] and virtual spectrophotometric [4] methods.

Table 5. The parameters used in the SPLITS model for each of the four sand samples, and the root-mean-square error (RMSE) between the actual and simulated reflectances.

View Table | View all tables in this article

One might suggest comparing SPLITS to existing models that have been used to represent sand. There are several problems with this approach. First, to the best of our knowledge, SPLITS is the first model proposed that provides both spectral and spatial (BRDF) responses for sand. Many existing models generate only the spatial response and use the spectral data as input. Furthermore, SPLITS uses sand characterization data as input. Candidates for comparison in the literature are generally designed to handle a wide variety of materials, and therefore use parameters not relating to sand. There is, therefore, no basis for comparison. Finally, even if two models yield similar responses, both models may be wrong. It is therefore best to compare directly with measured data.

These comparisons are made in two manners. First, we directly compare the output from the model to measured curves. Additionally, we demonstrate the effects caused by varying individual parameters to show that the resulting changes are in agreement with what is reported in the literature. We also compare the spatial response from the model to qualitative descriptions found in the literature. Additionally, images are presented in this section to illustrate the colours corresponding to the spectral responses obtained from the model.

5.1. Comparisons with Measured Data

Because the TEC database [75] contains a wide variety of sand samples and since this data covers the entire visible region of the spectrum, this data was chosen for comparison with the SPLITS model. Due to the lack of characterization data, however, precise quantitative comparisons are not feasible. Instead, the parameters were selected from within acceptable ranges (Table 2) as reported in the literature to attempt to obtain a good match between the measured and modeled curves. This allows us to demonstrate qualitative agreement with measured data, as well as to show that reflectance curves of actual sands are reproducible using the SPLITS model. Four samples were selected from the TEC database [75]: two dune sands, one from Australia (TEC #10019201), and one from Saudi Arabia (TEC #13_j9823), a magnetite rich beach sand from central Peru (TEC #10039240), and a sample from a dike outcrop in San Bernardino county, California (TEC #19au9815). In our experiments, we attempted to simulate the actual measurement conditions as accurately as possible according to the measurement setup outline provided by Rinker et al. [75].

5.1.1. Results of Comparisons

Figure 6 shows comparisons between the reflectance curves of these samples and the corresponding curves simulated using SPLITS. The reflectance is provided in the form of directional-hemispherical reflectance [54] with an incident angle of zero degrees. The corresponding parameters, along with root-mean-square (RMS) errors, are provided in Table 5. The low magnitude of the RMS errors, below 0.03, indicate a good spectral reconstruction even for remote sensing applications demanding high accuracy [6]. Computer generated images for each of the sand samples are depicted in Fig. 7. The sand grain pattern was extracted from a photograph and the spectral responses were provided by the SPLITS model. The degree of saturation, S, was also varied from S=0 to S=1 within each image to show the darkening effect simulated by the model.

 

Fig. 6. Comparisons between real and simulated sand. The solid line indicates the reflectance of the sand sample. The dashed line indicates the modeled reflectance. The parameters used and corresponding RMS errors are provided in Table 5. Top Left: TEC #10019201, Top Right: TEC #10039240, Bottom Left: TEC #13_j9823, Bottom Right: TEC #19au9815.

Download Full Size | PDF

 

Fig. 7. Computer generated images showing variation of sand colour with moisture as predicted by the model. The degree of saturation varies from S=0 at the top of each image to S=1 at the bottom. From left to right, the samples are TEC #10019201, TEC #10039240, TEC #13_j9823, and TEC #19au9815.

Download Full Size | PDF

 

Fig. 8. Qualitative behaviour of simulated sand. The solid line in each plot indicates the spectral response from the SPLITS model with ϑ_hg=ϑ_h+ϑ_g=0.01, r_hg=ϑ_h/(ϑ_h+ϑ_g)=0.00, ϑ_m=0.00, ζ₁=0.00, ζ₂=0.10, ζ₃=0.90, and µ′_p=µ′_m=µ′_c=1/3. Top Left: Varying ϑ_hg, Top Right: Varying r_hg, Bottom Left: Varying ϑ_m, Bottom Right: Varying the texture (ζ ₁,ζ₂,ζ₃). See Table 2 for symbol definitions.

Download Full Size | PDF

5.2. Qualitative Characteristics

To show that the SPLITS model qualitatively behaves like sand, we have also conducted simulations varying individual parameters. The variation in the spectral responses resulting from changes to these parameters are shown in Fig. 8. As expected [17], increasing the iron oxide concentration lowers the reflectance. Note also the shift in reflectance toward the red end of the spectrum when goethite is replaced with hematite (i.e., r_hg is increased), as confirmed in the literature [76, 22]. Also, as the particle size was decreased, the reflectance predicted by the model increased, in agreement with the literature [17, 28]. The resulting variation in colour is depicted in Fig. 9. Additionally, the degree of saturation varied between S=0 and S=1 for the four TEC sand samples. The darkening effect, reported in the literature [17, 48], is reflected in the SPLITS model. Simulated reflectance curves for two of these samples are given in Fig. 10.

A scattering simulation using SPLITS was also conducted to show the spatial distribution predicted by the model. This is shown in Fig. 11. While the predicted BRDF is diffuse, there is forward scattering and retro-reflection. This is also in agreement with the literature (Section 2.1).

 

Fig. 9. Computer generated images showing variation of sand colour as various parameters are changed. The image on the left (base image) corresponds to the solid lines in Fig. 8. The remaining images correspond to the spectral responses from the SPLITS model with the same parameters as in the base image except, from left to right, ϑ_hg=ϑ_h +ϑ_g=0.05, r_hg=ϑ_h/(ϑ_h+ϑ_g)=0.90, ϑ_m=0.30, ζ₃=1.0. See Table 2 for symbol definitions.

Download Full Size | PDF

 

Fig. 10. Spectral reflectance curves predicted by the SPLITS model for two of the TEC sand samples, varying the water content, expressed as the degree of saturation, S. Left: TEC #13_j9823, Right: TEC #19au9815.

Download Full Size | PDF

6. Conclusion

We have presented a new model that simulates spectral light transport for sand. The model applies a novel technique for simulating light interaction with particulate materials by generating these particles on the fly, rather than storing them explicitly. Although, in the context of this paper, we have applied this technique to the simulation of light transport with sand, the technique may potentially be applied to other heterogenous materials (organic and inorganic). Both the spectral and spatial responses (the measurement of appearance [55]) are represented using SPLITS. The model is controlled by parameters that relate to physical and mineralogical properties of sand, and the effects of these parameters on the model show good quantitative and qualitative agreement with observations reported in the literature.

While the results show good quantitative and qualitative agreement with measured data, they also indicate that there is still room for further improvement which is likely to depend on data availability. In fact, this paper also serves to point out the lack of and need for spectral BRDF measurements of natural surfaces along with the characterization data for those surfaces. Despite the efforts of numerous researchers [36, 39, 77], much work still needs to be done in this area. In most cases where the spectral reflectance or BRDF data is available, the corresponding characterization data is not, and vice versa. It is therefore difficult to verify quantitatively any models for these surfaces without making assumptions about the properties of the material in question. However, one can still validate the model qualitatively by demonstrating that modifying input parameters have the expected effect, and quantitatively by showing that it is possible to obtain matches between reflectances from real samples and from the model using plausible input parameters, as we have done.

 

Fig. 11. BRDF of simulated sand at 600nm. Three dimensional plots of the BRDF are shown on the top. Profiles along the principal plane are shown on the bottom (the dashed line indicates the incident direction). Left: Normal incidence, Middle: For light incident 30° from the normal, Right: For light incident 60° from the normal.

Download Full Size | PDF