In this paper, we present a new spectral light transport model for sand. The model employs a novel approach to simulate light interaction with particulate materials which yields both the spectral and spatial (bidirectional reflectance distribution function, or BRDF) responses of sand. Furthermore, the parameters specifying the model are based on the physical and mineralogical properties of sand. The model is evaluated quantitatively, through comparisons with measured data. Good spectral reconstructions were achieved for the reflectances of several real sand samples. The model was also evaluated qualitatively, and compares well with descriptions found in the literature. Its potential applications include, but are not limited to, applied optics, remote sensing and image synthesis.
© 2007 Optical Society of America
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 . 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  and virtual goniophotometric  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. , Brady , or Gerrard .
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) .
Soils are classified according to the size distribution of the mineral particles . 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 . In this work, we use the system developed by the United States Department of Agriculture (USDA) . 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).
2.1. Factors Affecting Light Transport
The reflectance of sand generally increases with wavelength . Additionally, sand reflection is non-Lambertian [15, 16]. In particular, sand exhibits retro-reflection , 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 .
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 .
2.1.1. Mineral Composition
As sand is composed primarily of weathered rock , 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  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 .
Iron oxide gives sand its distinctive hues. Indeed, the two most significant minerals that determine soil colour, hematite and goethite , are iron oxides. Goethite, also known as yellow ochre  or limonite , is one of the most common minerals found in soils . It colours soils yellow to brown . Hematite, or red ochre , imparts a red colour to soils and may mask the colour of goethite except when in small quantities . Hematite is usually found in tropical regions and is rare in temperate or cool climates . Iron oxides may be present as contaminants in the parent material . 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 .
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 . Smaller particles, and thus a higher density of scattering interfaces, result in higher reflectances [17, 18].
Grain shape may also affect scattering properties . 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 . In this paper, we use the sphericity measure proposed by Riley , 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 Di is the diameter of the largest inscribed circle and Dc is the diameter of the smallest circumscribed circle, as shown in Fig. 1.
In contrast, roundness can loosely be described as a measure of detail in the features on the grain surface . A higher roundness value indicates a smoother surface. Roundness is determined by comparing the particle to images on a roundness chart .
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 , Xie et al. , and Mishchenko et al. [34, 35].
3.1. General Models
Hapke  proposed an analytic radiative transfer model for particulate surfaces, which was then used to simulate reflectances from planetary surfaces . Emslie and Aronson  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 . Egan and Hilgeman  subsequently rectified this by simulating scattering, in addition to absorption, at these asperities. Oren and Nayar  generalized the widely used Lambertian model  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.  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  simulate light scattering in optically inhomogeneous particulate media using a geometric optics approach which combines ray tracing techniques and Monte Carlo methods. Recently, Peltoniemi  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
Various methods exist to achieve the darkening effect caused by water. Jensen et al.  use an extension of the Henyey-Greenstein phase function , 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 . It is also worth noting that it has been demonstrated by Mishchenko et al.  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  use an exponential function to empirically model reflectance in terms of moisture content. Neema et al.  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  apply Kubelka-Munk theory  to model the influence of iron oxides on soil colour. Torrent et al.  relate hematite and goethite concentrations to a function of the Munsell colour  of the soil, which Torrent et al. call the redness rating. Okin and Painter  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 ) responses are simulated. These two components together constitute the measurement of appearance . 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” .
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.
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 interference1, 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 . 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 .
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 .
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 . This may be related to the degree of saturation by S=w/P.
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 , 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 . 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
where X is normally distributed with zero mean and standard deviation
This standard deviation was chosen so that R<n′·n≤1 for 95% of the facets . 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 .
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  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 . Recently, it was applied to assess the effect of grain size on the reflectance of sandy materials . 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 .
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
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
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  then yields
The volume of the coating may be estimated as ASh, where AS is the surface area of the particle, since we are assuming that the h≪s. Therefore,
where V is the volume of the particle. Dividing the numerator and denominator by V yields
where AV (s,Ψ)=AS(s,Ψ)/V(s,Ψ) is the surface area to volume ratio of the particle. Noting that
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 
For simplicity, we use the mean sphericity, , in the above (i.e., AV (1, )) 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 , 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. . 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 dg and σg for each soil separate. Defining
where ag=logdg and bg=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 . Shirazi et al.  provide a table for determining dg 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 , 2 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  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.
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/Kj, where Kj is the geometric attenuation coefficient for particles of type j. This distance is given by
where ξ is uniformly distributed on (0,1) .
For prolate spheroids with normally distributed sphericity and log-normally distributed size, Kj reduces  to
The particle size distribution, fg(s), is given by Eq. (19), AV (1,Ψ) is given by Eq. (14), and Φ′(x̄,σ2)(x) is the probability density function for the normal distribution with mean x̄ and variance σ2 .
The particle size, s, is then chosen according to the probability density function
and the sphericity, Ψ, is chosen according to
where C 1 and C 2 are the constants that ensure that fs and fΨ (respectively) integrate to one over their domain . Note that Eq. (25) does not match the distribution provided by Shirazi et al.  shown in the integrand of Eq. (20). This follows from the distinction between the particle size distribution by mass, provided by Shirazi et al.  (or by volume as we are interpreting it), and the distribution of particle sizes struck by rays travelling randomly through the medium .
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-1), where ξ1 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 is the eccentricity of spheroid . 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πξ 2 (ξ2 being another canonical random variable) and . 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  (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 . The ray is then reflected in n′ with a probability once again determined using the Fresnel equations3 , considering the media on either side of the interface. Otherwise, the ray is refracted according to Snell’s Law for absorbing media . 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) .
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 . The light is transmitted a distance d with probability T=exp(-αd), where α is the absorption coefficient, given by
. 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.
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 . 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 , 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 .
To evaluate the SPLITS model objectively, it is necessary to perform comparisons with measured data directly. Thus, comparisons were made between measured data  and results obtained from the model using virtual goniophotometric  and virtual spectrophotometric  methods.
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  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 : two dune sands, one from Australia (TEC #10019201), and one from Saudi Arabia (TEC #13j9823), 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. .
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  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 . 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.
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 , 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., rhg 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).
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 ) 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.
The authors would like to thank the anonymous reviewers for their helpful comments. The work presented in this paper was partially supported by the Natural Sciences and Engineering Research Council of Canada (NSERC grant 238337) and the Canada Foundation for Innovation (CFI grant 33418).
|1||It had been demonstrated that the Monte Carlo simulation of the optical path of light propagating through n scattering events in an inhomogeneous medium is equivalent to the calculation of the n th order ladder diagram, and the Monte Carlo approach can be extended to the computation of time correlation functions and coherent backscattering .|
|2||Recall that the extended boundary represents the interface between the pore space (air or water) and the ambient medium. It does not represent the rough sand surface. This arises from the random positioning of the simulated particles.|
|3||Note that the Fresnel equations are applied to the microfacets, not the overall particle surface. Since the microfacets themselves are planar surfaces, the Fresnel equations may be applied.|
References and links
1. S. Prahl, “Light Transport in Tissue,” Ph.D. thesis, University of Texas at Austin (1988).
2. S. Prahl, M. Keijzer, S. Jacques, and A. Welch, “A Monte Carlo Model of Light Propagation in Tissue,” SPIE Institute Series 5, 102–111 (1989).
4. G. Baranoski, J. Rokne, and G. Xu, “Virtual Spectrophotometric Measurements for Biologically and Physically Based Rendering,” The Visual Computer 17, 506–518 (2001). [CrossRef]
5. A. Krishnaswamy, G. Baranoski, and J. G. Rokne, “Improving the Reliability/Cost Ratio of goniophotometric measurement,” 9, 31–51 (2004).
6. S. Jacquemoud, S. Ustin, J. Verdebout, G. Schmuck, G. Andreoli, and B. Hosgood, “Estimating Leaf Biochemistry Using PROSPECT Leaf Optical Properties Model,” Remote Sensing of Environment 56, 194–202 (1996). [CrossRef]
7. R. Shuchman and D. Rea, “Determination of Beach Sand Parameters Using Remotely Sensed Aircraft Reflectance Data,” Remote Sensing of Environment 11, 295–310 (1981). [CrossRef]
8. R. Morris and D. Golden, “Goldenrod Pigments and the Occurrence of Hematite and Possibly Goethite in the Olympus-Amazonis Region of Mars,” Icarus 134, 1–10 (1998). [CrossRef]
9. R. Singer, “Spectral Evidence for the Mineralogy of High-Albedo Soils and Dust on Mars,” Journal of Geophysical Research 87, 10,159–10,168 (1982). [CrossRef]
10. F. Pettijohn, P. Potter, and R. Siever, Sand and Sandstone, 2nd ed. (Springer-Verlag, New York, NY, 1987). [CrossRef]
11. N. Brady, The Nature and Properties of Soils, 8th ed. (Macmillan Publishing Co., Inc., New York, NY, 1974).
12. J. Gerrard, Fundamentals of Soils (Routledge, New York, NY, 2000).
13. Soil Science Division Staff, Soil Survey Manual (Soil Conservation Service, 1993). United States Department of Agriculture Handbook 18.
14. K. Coulson and D. Reynolds, “The Spectral Reflectance of Natural Surfaces,” Journal of Applied Meteorology 10, 1285–1295 (1971). [CrossRef]
15. J. Norman, J. Welles, and E. Walter, “Contrasts Among Bidirectional Reflectance of Leaves, Canopies, and Soils,” IEEE Transactions on Geoscience and Remote Sensing GE-23, 659–667 (1985). [CrossRef]
16. K. Coulson, G. Bouricius, and E. Gray, “Optical Reflection Properties of Natural Surfaces,” Journal of Geophysical Research 70, 4601–4611 (1965). [CrossRef]
17. M. Baumgardner, L. Silva, L. Biehl, and E. Stoner, “Reflectance Properties of Soils,” Advances in Agronomy 38, 1–43 (1985). [CrossRef]
18. D. Leu, “Visible and Near-Infrared Reflectance of Beach Sands: A Study on the Spectral Reflectance/Grain Size Relationship,” Remote Sensing of Environment 6, 169–182 (1977). [CrossRef]
19. G. Hunt and J. Salisbury, “Visible and Near-Infrared Spectra of Minerals and Rocks: I. Silicate Minerals,” Modern Geology 1, 283–300 (1970).
20. P. Farrant, Color in Nature: A Visual and Scientific Exploration (Blandford Press, 1999).
21. A. Mottana, R. Crespi, and G. Liborio, Simon and Schuster’s Guide to Rocks and Minerals (Simon and Schuster, Inc., New York, NY, 1978).
22. J. Torrent, U. Schwertmann, H. Fechter, and F. Alferez, “Quantitative Relationships Between Soil Color and Hematite Content,” Soil Science 136, 354–358 (1983). [CrossRef]
23. R. Cornell and U. Schwertmann, The Iron Oxides, 2nd ed. (Wiley-VCH GmbH & Co. KGaA, Weinheim, Germany, 2003). [CrossRef]
24. H. Wopfner and C. Twindale, “Formation and Age of Desert Dunes in the Lake Eyre Depocentres in Central Australia,” Geologische Rundschau 77, 815–834 (1988). [CrossRef]
25. G. Hunt, J. Salisbury, and C. Lenhoff, “Visible and Near-Infrared Spectra of Minerals and Rocks: III. Oxides and Hydroxides,” Modern Geology 2, 195–205 (1971).
26. S. Twomey, C. Bohren, and J. Mergenthaler, “Reflectance and Albedo Differences Between Wet and Dry Surfaces,” Applied Optics 25, 57–84 (1986). [CrossRef]
27. M. Kühl and B. Jørgensen, “The Light Field of Microbenthic Communities: Radiance Distribution and Microscale Optics of Sandy Coastal Sediments,” Limnology and Oceanography 39, 1368–1398 (1994). [CrossRef]
29. H. Wadell, “Volume, Shape, and Roundness of Rock Particles,” Journal of Geology 40, 443–451 (1932). [CrossRef]
30. N. Riley, “Projection Sphericity,” Journal of Sedimentary Petrology 11, 94–95 (1941).
31. W. Krumbein, “Measurement and Geological Significance of Shape and Roundness of Sedimentary Particles,” Journal of Sedimentary Petrology 11, 64–72 (1941).
32. H. Zhang and K. Voss, “Comparisons of Bidirectional Reflectance Distribution Function Measurements on Prepared Particulate Surfaces and Radiative-Transfer Models,” Applied Optics 44, 597–610 (2005). [CrossRef] [PubMed]
33. Y. Xie, P. Yang, B.-C. Gao, G. Kattawar, and M. Mishchenko, “Effect of Ice Crystal Shape and Effective Size on Snow Bidirectional Reflectance,” Journal of Quantitative Spectroscopy and Radiative Transfer 100, 457–469 (2006). [CrossRef]
34. M. Mishchenko, L. Travis, and A. Lacis, Multiple Scattering of Light by Particles: Radiative Transfer and Coherent Backscattering (Cambridge University Press, Cambridge, 2006).
36. B. Hapke, “Bidirectional Reflectance Spectroscopy. 1. Theory,” Journal of Geophysical Research 86, 3039–3054 (1981). [CrossRef]
37. B. Hapke and E. Wells, “Bidirectional Reflectance Spectroscopy. 2. Experiments and Observations,” Journal of Geophysical Research 86, 3055–3054 (1981). [CrossRef]
40. M. Oren and S. Nayar, “Generalization of Lambert’s Reflectance Model,” in Computer Graphics Proceedings, Annual Conference Series, pp. 239–246 (1994).
41. L. Wolff, “Diffuse-Reflectance Model for Smooth Dielectric Surfaces,” Journal of the Optical Society of America A (Optics, Image Science, and Vision) 11, 2956–2968 (1994). [CrossRef]
42. M. Mishchenko, J. Dlugach, E. Yanovitskij, and N. Zakharova, “Bidirectional reflectance of flat, optically thick particulate layers: An efficient radiative transfer solution and applications to snow and soil surfaces,” Journal of Quantitative Spectroscopy and Radiative Transfer 63, 409–432 (1999). [CrossRef]
43. D. Stankevich and Y. Shkuratov, “Monte Carlo Ray-Tracing Simulation of Light Scattering in Particulate Media with Optically Contrast Structure,” Journal of Quantitative Spectroscopy and Radiative Transfer 87, 289–296 (2004). [CrossRef]
44. J. Peltoniemi, “Spectropolarised Ray-Tracing Simulations in Densely Packed Particulate Medium,” Journal of Quantitative Spectroscopy and Radiative Transfer (2007). In Press, accepted, URL http://dx.doi.org/10.1016/j.jqsrt.2007.05.009. [CrossRef]
45. H. Jensen, J. Legakis, and J. Dorsey, “Rendering ofWet Materials,” in Proceedings of the Eurographics Workshop on Rendering, pp. 273–282 (1999).
46. L. Henyey and J. Greenstein, “Diffuse Radiation in the Galaxy,” Astrophysical Journal 93, 70–83 (1941). [CrossRef]
47. Z. Li, A. Fung, S. Tjuatja, D. Gibbs, C. Betty, and J. Irons, “A Modeling Study of Backscattering from Soil Surfaces,” IEEE Transactions on Geoscience and Remote Sensing 34, 264–271 (1996). [CrossRef]
48. D. Lobell and G. Asner, “Moisture Effects on Soil Reflectance,” Soil Science Society of America Journal 66, 722–727 (2002). [CrossRef]
49. D. Neema, A. Shah, and A. Patel, “A Statistical Optical Model for Light Reflection and Penetration Through Sand,” International Journal of Remote Sensing 8, 1209–1217 (1987). [CrossRef]
50. V. Barron and L. Montealegre, “Iron Oxides and Color of Triassic Sediments: Application of the Kubelka-Munk Theory,” American Journal of Science 286, 792–802 (1986). [CrossRef]
51. P. Kubelka and F. Munk, “Ein Beitrag zur Optik der Farbanstriche (An Article on Optics of Paint Layers),” Zeitschrift für Technische Physik 12, 593–601 (1931).
52. D. Nickerson, “History of the Munsell Color System and its Scientific Application,” Journal of the Optical Society of America 30, 575–645 (1940). [CrossRef]
53. G. Okin and T. Painter, “Effect of Grain Size on Remotely Sensed Spectral Reflectance of Sandy Desert Surfaces,” Remote Sensing of Environment 89, 272–280 (2004). [CrossRef]
54. F. Nicodemus, J. Richmond, J. Hsia, I. Ginsberg, and T. Limperis, Geometrical Considerations and Nomenclature for Reflectance (National Bureau of Standards, United States Department of Commerce, 1977).
55. R. Hunter and R. Harold, The Measurement of Appearance, 2nd ed. (JohnWiley and Sons, New York, NY, 1987).
56. G. Ward, “Measuring and Modeling Anisotropic Reflection,” Computer Graphics 26, 262–272 (1992). [CrossRef]
57. V. Kuz’min and I. Meglinski, “Numerical Simulation of Coherent Effects under Conditions of Multiple Scattering,” Optics & Spectroscopy 97, 100–106 (2004). [CrossRef]
58. F. Pedrotti and L. Pedrotti, Introduction to Optics, 2nd ed. (Prentice Hall, Upper Saddle River, NJ, 1993).
59. T. Nousiainen, K. Muinonen, and P. Räisänen, “Scattering of Light by Large Saharan Dust Particles in a Modified Ray Optics Approximation,” Journal of Geophysical Research 108, AAC 12–1–17 (2003). [CrossRef]
61. W. Tropf, M. Thomas, and T. Harris, “Properties of Crystals and Glasses,” in Handbook of Optics, M. Bass, E. Van Stryland, D. Williams, and W. Wolfe, eds., vol. 2, 2nd ed., chap. 33 (McGraw-Hill, 1995).
62. I. Sokolik and O. Toon, “Incorporation of Mineralogical Composition into Models of the Radiative Properties of Mineral Aerosol from UV to IR Wavelengths,” Journal of Geophysical Research 104, 9423–9444 (1999). [CrossRef]
63. W. Egan and T. Hilgeman, Optical Properties of Inhomogeneous Materials: Applications to Geology, Astronomy, Chemistry, and Engineering (Academic Press, New York, NY, 1979). [PubMed]
64. A. Schlegel, S. Alvarado, and P. Wachter, “Optical Properties of Magnetite (Fe3O4),” Journal of Physics C: Solid State Physics 12, 1157–1164 (1979). [CrossRef]
65. M. Vepraskas and D. Cassel, “Sphericity and Roundness of Sand in Coastal Plain Soils and Relationships with Soil Physical Properties,” Soil Science Society of America Journal 51, 1108–1112 (1987). [CrossRef]
66. B. Kimmel, “SPLITS: A Spectral Light Transport Model for Sand,” Master’s thesis, School of Computer Science, University of Waterloo (2005).
67. S. Ross, A First Course in Probability, 5th ed. (Prentice Hall, Upper Saddle River, NJ, 1998).
68. C. Bohren and D. Huffman, Absorption and Scattering of Light by Small Particles (John Wiley and Sons, New York, NY, 1983).
69. J. Maxwell Garnett, “Colours in Metal Glasses and in Metallic Films,” Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character 203, 385–420 (1904). [CrossRef]
70. M. Shirazi, L. Boersma, and J. Hart, “A Unifying Quantitative Analysis of Soil Texture: Improvement of Precision and Extention of Scale,” Soil Science Society of America Journal 52, 181–190 (1988). [CrossRef]
71. J. Snyder, Map Projections: A Working Manual (United States Government Printing Office, Washington, 1987). U.S. Geological Survey Professional Paper 1395.
72. B. Hapke, “Bidirectional Reflectance Spectroscopy. 4. The Extinction Coefficient and the Opposition Effect,” Icarus 67, 264–280 (1986). [CrossRef]
73. M. Born and E. Wolf, Principles of Optics, 6th ed. (Pergamon Press, Oxford, England, 1980).
74. C. Gribble and A. Hall, Optical Mineralogy: Principles and Practice (Chapman & Hall, New York, NY, 1993).
75. J. Rinker, C. Breed, J. McCauley, and P. Corl, “Remote Sensing Field Guide — Desert,” Tech. rep. , U.S. Army Topographic Engineering Center, Fort Belvoir, VA (1991).
76. T. Cudahy and E. Ramanaidou, “Measurement of the Hematite:Goethite Ratio Using Field Visible and Near-Infrared Reflectance Spectrometry in Channel Iron Deposits, Western Australia,” Australian Journal of Earth Sciences 44, 411–420 (1997). [CrossRef]
77. J. Cierniewski, “A Model for Soil Surface Roughness Influence on the Spectral Response of Bare Soils in the Visible and Near-Infrared Range,” Remote Sensing of Environment 23, 92–115 (1987). [CrossRef]