Background

DIRSIG4 took an open-ended approach to optical properties — individual models were independently implemented and were only required to respond to basic model queries. This approach was particularly well suited for an academic environment since new models could be brought in directly and easily, often by simply translating equations in papers into code. However, that same open-endedness meant that the level of initialization, general computational complexity, and the potential ability to optimize, varied widely depending on the user’s choice of model. More complex optical property models also needed to be appropriately paired with radiometry solvers that could use the properties as effectively as possible — often requiring fine-tuning of many detailed runtime parameters to achieve both radiometric accuracy and optimal speed.

Our approach in DIRSIG5 was to unify the diversity of optical property model options into a three-prong system. We would provide optimized implementations of "ideal" optical properties (perfectly diffuse, "Lambertian" properties as well as singularly valued or "delta" functions) and a single data-based model to handle everything else. In the development process we added a fourth option to provide an intermediate step above fully diffuse that allows the definition of a specular lobe without having to go all the way to a generic data model for some commonly used, simple models.

As with the models themselves, our approach to optical property radiometry was to replace the various parametrized solvers in DIRSIG4 with a single, unified approach based on Monte Carlo integration. This change means that the quality of radiometry is driven by a single, global parameter (the criteria used to converge to a solution) rather than the myriad solver-specific options that were in place previously.

The Data Model

The data model is a catch-all for any optical properties that are not handled directly with the aforementioned simple analytical models and the way to bring generic optical property data directly into DIRSIG5. The first component of the model is a high spectral resolution integral of the function (i.e. the directional hemispherical reflectance, or DHR), which is stored and handled in a way that is similar to the Lambertian model. The second component is a shaping function that distributes energy into the hemisphere or sphere, i.e. it defines the PDF (probability density function) of the property. It is this function that drives the critical sampling function during a simulation and, in combination with the integral, the direct contribution of sources such as the sun.

A number of different approaches to the PDF model were attempted, including projecting onto basis function representations such as Zernike polynomials and spherical harmonics. However, accurate representation of non-trivial properties required a significant number of coefficients, especially when oscillating on tight lobes that don’t align with the underlying basis functions.

In the end, we decided to take a data-based approach based on a partition of the hemisphere (sphere). The partition is based on the Hierarchical Equal Area isoLatitude Pixelization, or HEALPix, subdivision of the sphere developed by NASA/JPL (https://healpix.jpl.nasa.gov). Each cell of the base, 12 diamond-shaped cell partition over the sphere is successively split into four sub-cells with each depth level step in the hierarchy. The number of cells for a given partition is equal to \( 12 \cdot 4^d\), where d is the depth of the partition. Examples of the parition for depths 3, 4, and 5 are shown below. The corresponding number of cells at each depth are 768, 3072, and 12288, respectively.

sh
Figure 1. Example HEALPix sphere partitions for the given depth.

While there are many ways to partition the sphere, the HEALPix partition is notable for maintaining equal area cells and having a consistent, compact shape simultaneously (it is not possible to subdivide the sphere into cells that have equal area and equally spaced centers, so there is always a tradeoff). The former is essential for sampling (we don’t want to do on-the-fly corrections for varying solid angles), and the latter for appropriately localizing a cell on the sphere. While public source code for the official HEALPix partition is available from JPL, its (GPL) license is not compatible with DIRSIG and our implementation has be developed independently based on:

Roukema, B. F., & Lew, B. (2004). A Solution to the isolatitude, equi-area, hierarchical pixel-coordinate system. arXiv preprint astro-ph/0409533.

As such, indexing and alignment may vary between the DIRSIG5 implementation and the HEALPix library.

Our usage of the HEALPix partition takes advantage of the hierarchical nature of the subdivision, treating the representation as twelve quad-tree data structures (one quad-tree for each of the original twelve cells). This allows for quick traversal of the data (the maximum number of steps equal to the full depth) and we do an initial data compression based on comparison of child-nodes (i.e. approximately constant valued child cells can be collapsed into the parent and less data stored). Finally, we can use the diamond-shape structure itself as the basis for an interpolation scheme across the surface of the sphere to overcome quantization.

harmonics
Figure 2. Result of randomly sampling an SQT-stored "BRDF" (right) constructed with data from a (non-physical) set of spherical harmonics generated on a spherical grid (left). Apparent "holes" in the queried data are areas which were not sampled (i.e. there were no queries at those points).

Terminology

BRDF

The Bidirectional Reflectance Distribution Function (or BRDF) is the most common optical property defined by a user for DIRSIG. It describes how incident planar irradiance is scattered off of a surface per unit area as an exitant radiance.

\begin{equation*} \rho_{r,\lambda}(\nu_e,\nu_i)=\frac{dL_\lambda(\omega_e)}{dE_\lambda(\omega_i)}=\ \frac{dL_\lambda(\omega_e)}{L_\lambda(\omega_i)(n\cdot \nu_i)d\omega_i}, \end{equation*}

where the incident planar irradiance has been expanded to show that it is equal to the projected radiance from the differential solid angle around the incident direction (the projection is a cosine effect equal to the dot product between the normal vector and the incident direction and is a consequence of maintaining a constant area).

DHR

The Directional-Hemispherical Reflectance (or DHR) is constructed from the BRDF by integrating over all possible directions in the hemisphere and accounting for the area projections:

\begin{equation*} \bar\rho_{r,\lambda}(\nu_e)=\int_{\Omega_i} \rho_{r,\lambda}(\nu_e,\nu_i) (n\cdot \nu_i)d\omega_i, \end{equation*}

where the maximum value of the DHR is unit valued in order to conserve energy (in contrast, the maximum value of the BRDF can be much larger than one). While the DHR is useful on its own to understand the behavior of a surface (especially when compared to the transmittance equivalent), we can also think of it as the magnitude of the cosine-weighted BRDF function for a given exitant vector.

Decomposition of the Rendering Equation

When we encounter a surface in DIRSIG we are usually trying to calculate the radiance leaving that surface into a particular exitant direction. Ignoring self-emission, the so-called "rendering equation" is the integral form of the BRDF definition and gives the surface-leaving radiance as:

\begin{equation*} L_{\lambda}(\nu_e)=\int_{\Omega_i}\rho_{r,\lambda}(\nu_e,\nu_i)(n\cdot \nu_i)L_{\lambda}(\nu_i) d\omega_i. \end{equation*}

Using the concept of the spectral DHR as the magnitude of the cosine-weighted BRDF, we can rewrite this function as

\begin{equation*} L_{\lambda}(\nu_e)=\bar\rho_{r,\lambda}\int_{\Omega_i}P_{\bar\rho,\lambda}(\nu_e,\nu_i)L_{\lambda}(\nu_i) d\omega_i, \end{equation*}

or, in other words, the DHR-weighted integral over the probability density function (PDF) of the cosine-weighted BRDF.

The integral portion is approximately equivalent to the expected value of the function when samples are drawn from the PDF (i.e. we’ve constructed a Monte Carlo integral). Using this relation allows us to remove the explicit integral entirely,

\begin{equation*} L_{\lambda}(\nu_e)\approx\bar\rho_{r,\lambda}\mathbb{E}\left\{L_{\lambda}(\hat\nu_i)\right\}, \end{equation*}

or the average value of the radiance obtained by sampling the incident direction from the PDF, then scaled by the DHR. The quality of the approximation then goes up with the number of samples that are used.

This equation, and others like it, are the foundation of tbe path-based radiometry in DIRSIG5. The integral is solved indirectly and over an aggregation of samples for a particular surface and it is not necessary to use the BRDF directly (though direct queries are necessary for incorporating direct sources, such as the sun or moon).

Structure of the SQT

In order to provide the functionality described above (BRDF and DHR queries plus sampling), the structure of the SQT is as optimized as possible. First, we recognize that incorporating the projected area is essential for maximizing the contribution of each sample. So rather than store the BRDF directly, we store the cosine weighted BRDF for a given exitant angle.

Second, the primary role of the SQT is to provide a quick route to sampling the cosine weighted BRDF. We pre-normalize the hemisphere so we are left with the PDF of the DHR. This allows us to traverse the 12 quad-trees with a single random number draw from [0,1]. The scalar magnitude of the DHR is stored separately and can be recombined when querying the (cosine-weighted) BRDF.

This separation of the shape from the magnitude suggest our approach to handling spectral variation. Each SQT is only valid for a single wavelength (and a single source/view vector) and it would be very costly to rely on an independent SQT for every posssible wavelength in a simulation. Instead we recognize that in most instances, the "shape" (i.e. the PDF) of the reflectance function does not vary as quickly as the magnitude (i.e. the magnitude). Therefore, we can define the shape at a coarser spectral sampling rate than the magnitude and simply interpolate the shape between spectral samples. This allows us to bring in more finely sampled DHR data than we might have for full bidirectional data.

RAW file format

The "RAW" optical format is intended to hold sample optical property data for one or more exitant angles, a set of wavelength samples, and arbitrary scattered directions. Each set of optical property data is stored in a single file containing:

  • a header with an 8 character signature and any metainfo

  • a list of wavelengths at which the property was sampled

  • a list of exitance angles at which the property was sampled (only for bidirectional properties); if the property is anisotropic, then there are both declination and azimuth angles

  • spectral data for each sample (in wavelength list order), preceded by the exitant angle index (if bidrectional) and the vector corresponding to the sample direction

The raw data can be irregularly spaced across the hemisphere or sphere. The exitant angles and wavelength samples can’t change though and data must be provided for each spectral point.

RAW header

The header of the RAW file is a human-readable ASCII line of 256 characters. This allows anyone to quickly view the contents of a file whether or not the data is stored in an ASCII or binary format.

The first 8 characters of the header form the signature and are used to indicate the (abstracted) type of optical property being described as well as the version and format of the RAW file, and any additional author/descriptive information. All 256 characters (including whitespace) must be used for the binary format, but can be truncated in the ASCII version. The layout is as follows:

R

A

W

B/U/A

H/S

1

0

A/B

metainfo

0

1

2

3

4

5

6

7

8-255

The meaning of each location is explained below:

[0-2]

"RAW" indicates the basic type of file that follows

[3]

if "B": the property is bidirectional (e.g. a BRDF)
if "U": the property is unidirectional (e.g. a SPF)
if "A": the property is bidirectional and also anisotropic
(e.g. a BRDF that varies with exitance azimuth angle)

[4]

if "H": the data is hemispherical only (e.g. BRDF)+ if "S": the data is spherical (e.g. BSDF)

[5]

version (major)

[6]

version (minor)

[7]

if "A": data is stored in ASCII format
if "B": data is stored in a binary format

[8-255]

arbitrary information about the data file
(e.g. author, date, description)

RAW sub-header: wavelengths

All types of the RAW format define the wavelengths at which the subsequent data was measured. The sub-header consists of two parts. First the number of wavelengths are given, then the actual samples.

Nλ

λ1

…​

λN

uint32

float32

…​

float32

Numbers are space delimited and followed by a linebreak in the ASCII version (no formatting in the binary version).

All wavelengths are given in units of Microns (μm). They can be in non-sequential order and arbitrarily spaced. If the same wavelength is given multiple times, the behavior is undefined.

RAW sub-header: exitance declination angles (bidirectional properties only)

Data in the RAW format defining a bidrectional property (e.g. BRDF) requires a sub-header defining the exitance declination angles that drive the measurements (i.e. the angle from nadir at which the surface is observed). Note that if the property being described is isotropic then a single angle is sufficient to define the viewer (the azimuth angle section below is not used). The source is assumed to be in the X-Z plane and the exitant angle is always measured towards +X in units of radians.

As with the wavlengths sub-header, the exitance declination angles consist of two components, the number of angles and then the list.

Nθ

θ1

…​

θN

uint32

float32

…​

float32

Values should always be in the range [0:π] (though BRDF values would only be to π/2). Numbers are space delimited and followed by a linebreak in the ASCII version (no formatting in the binary version).

RAW sub-header: exitance azimuth angles (bidirectional, anisotropic properties only)

Data in the RAW format defining an anisotropic bidrectional property (e.g. an anisotropic BRDF) requires a sub-header defining the exitance azimuthal angles that drive the data as measured CCW from an arbitrary reference direction (an azimuth of zero). Regardless of azimuth, the viewer is still assumed to be in the X-Z (principal) plane towards +X for purposes of defining source vectors (internally, the correct data set will be rotated as needed). This version of the RAW file assumes symmetry across the principal plane (i.e. the data for +φ is the same as -φ)

The exitance azimuth angles consist of two components, the number of angles and then the list. The first angle in the list must be zero and the maximum is $\pi$.

Nφ

φ1

…​

φN

uint32

float32

…​

float32

Values should always be in the range [0:π] (anisotropic properties are still assumed to be symmetric across the principal plane). Numbers are space delimited and followed by a linebreak in the ASCII version (no formatting in the binary version).

All isotropic functions should skip this input (the implicit azimuth is zero).

RAW Spectral Data

The remainder of the data is provided as spectral BRDF data at exitant angle indexes (0-based indexing corresponding to the list of exitant angles above) The exitant declination index is provided first, then the azmiuth index (only if the property is anisotropic), then the 3-element vector (unit length) corresponding to the BRDF measurement location. Finally, one value is provided for each of the spectral samples listed in the header. Spectral values should be non-negative, finite numbers and should not include the cosine due to the projected area.

indexθ

indexφ**

vx

vy

vz

value1

…​

value

uint32

uint32**

float32

float32

float32

float32

…​

float32

** should be omitted unless the property is anisotropic

SQT file format

The "SQT" optical format is intended to hold preprocessed optical property data that has been sorted into a Spherical Quad-Tree. Besides the human-readable header, it is exclusively machine generated, binary, may be compressed internally. It is used to quickly load optical property data into a simulation and should only be inspected with appropriate tools (though the header is human readable). It consists of:

  • a 1024 character ASCII header with an 8 character signature and any metainfo

  • binary, compressed data for the SQT (internal format, version dependent)

SQT header

The header of the SQT file is a human-readable ASCII line of 1024 characters. This allows anyone to quickly view the contents of a file.

The first 8 characters of the header form the signature and are used to indicate the (abstracted) type of optical property being described the version and format of the SQT file, and any additional author/descriptive information. Since it is a fixed binary format, all 1024 characters are stored, regardless of whether they are all used.

S

Q

T

B/U/A

H/S

1

0

R/S

metainfo

0

1

2

3

4

5

6

7

8-255

The meaning of each location is explained below:

[0-2]

"SQT" indicates the basic type of file that follows

[3]

if "B": the property is bidirectional (e.g. a BRDF)
if "U": the property is unidirectional (e.g. a SPF)
if "A": the property is bidirectional and also anisotropic
(e.g. a BRDF that varies with exitance azimuth angle)

[4]

if "H": the data is hemispherical only (e.g. BRDF)+ if "S": the data is spherical (e.g. BSDF)

[5]

version (major)

[6]

version (minor)

[7]

if "R": data was generated by converting a RAW file
if "D": data was automatically generated by the DIRSIG 4 scene compilation tool (scene2hdf)

[8-1023]

arbitrary information about the data file
(e.g. author/label, date, description)