This article explains how to write a DLL for a user defined scattering function. An example of a Gaussian X-Y probability function is given.
Authored By Sanjay Gangadhara
Introduction
Light scattered from surfaces can contribute significantly to the image quality and overall performance of an optical system. Oftentimes, it is a desired effect which is why OpticStudio provides Lambertian, Gaussian, and ABg surface scattering models. However, there may be cases where it is necessary to model objects whose surface scattering properties do not follow any of the built-in distributions. In such cases, user-defined scattering models may be employed via a DLL.
This article will provide some of the theory behind surface scattering, and will show an example of using a DLL to create a user-defined scatter distribution.
DLL help
A DLL is a Windows Dynamic Link Library (DLL), as described in the Knowledgebase article “How to compile a User-Defined DLL”. Please read through and understand that article before proceeding with this one. If you have the need for a user defined scattering model but do not wish to write the software yourself, we have a list of OpticStudio consultants on our website.
There are a few examples of user-defined surface scattering models which are provided with your OpticStudio installation. Both the C source code and corresponding DLL file for each model are located in the directory Documents\Zemax\DLL\SurfaceScatter (assuming Opticstudio 17+). The majority of each code is similar, containing “boilerplate” programming used to link between the DLL and OpticStudio. The difference between the models is in the specification of the direction cosines for the scattered ray. This article describes how to determine those values for any analytic, integrable scattering distribution.
Generating the statistical description of a scattering distribution described by probability
The direction cosines of any scattered ray are determined using a Monte Carlo technique. Two numbers (RN1 and RN2) are drawn from a random number generator designed to produce values between 0 and 1 with uniform probability. These numbers then define two direction cosines for the scattered ray. The third direction cosine is determined from the constraint that the magnitude of the scattered ray vector is 1.
The question is then how the direction cosines are related to the random numbers. This relationship is determined from integration of the scattering distribution. Suppose that the scattering model is given in terms of a probability, e.g. say that P(θ,ϕ) describes the probability of the ray scattering into a range of angles dθ (centered about θ) and dϕ (centered about ϕ):
where θ is the polar angle with respect to the surface normal n and ϕ is the azimuthal angle, as shown above. The necessary integral, performed in polar coordinates, is then:
This indefinite integral (which is still a function of θ and ϕ) is then normalized by the Total Integrated Scatter (TIS). The TIS is calculated using the exact same formula given above, but it is a definite integral, i.e. the integration is performed over the full limits of θ and ϕ (i.e. the integration goes from 0 to π/2 in θ and from 0 to 2π in ϕ). The normalized integral (IN = I/TIS) is then a function of θ and ϕ which – when taking into account the constant of integration that arises from the indefinite integration P(θ,ϕ) – is constrained to be between 0 and 1.
For the case when P(θ,ϕ) is separable in θ and ϕ – i.e. P(θ,ϕ) = F(θ)*G(ϕ) – IN can be separated into two terms:
Values for the scattered angles are determined by setting each of the two random numbers equal to Iθ and Iϕ, respectively, and then inverting the equations to solve for θ and ϕ:
If P(θ,ϕ) is not separable, then an intermediate step is necessary, but the basic idea is the same. Iθ(θ,ϕ) is first calculated by integrating over ϕ:
This indefinite integral is normalized to the definite integral (= Iθ(θ)) calculated by performing the same integration over the full limits of ϕ (0 to 2p). I and TIS are then calculated by performing the indefinite and definite integration of Iθ(θ), respectively:
The value for the polar angle is determined by setting the first random number equal to IN and inverting the equation to solve for θ:
The value for the azimuthal angle is then determined by setting the second random number equal to Iθ(θ,ϕ) for the known value of θ, and inverting to solve for ϕ:
For this case when P(θ,ϕ) is not separable, the order of integration is in fact not important. We could have chosen to instead integrate P(θ,ϕ) over θ first, and then over ϕ. The technique for relating the random numbers RN1 and RN2 to θ and ϕ would have been identical:
Integration may therefore be performed in whichever order is most convenient for the given probability function.
Generating the statistical description of a scattering distribution described by BSDF
If the scattering model is provided in terms of the Bi-Directional Scatter Distribution Function (BSDF), then the integrals necessary to relate the random numbers to the scatter angles include an extra cos(θ) factor:
The cosine factor arises from the fact that BSDF is defined in terms of the irradiance falling on the surface, and thus requires use of the projected surface area (Ap = A*cos(θ)). Otherwise, the procedure for determining the scatter angles θ and ϕ is identical to the procedure described in the previous section.
Let’s look at a simple example – the case of a Lambertian scattering source. For such a source, the BSDF is constant for all angles, equal to 1/π. I is then given by:
The total integrated scatter is calculated by setting the limits of integration from 0 to π/2 in θ and from 0 to 2π in ϕ:
Therefore I is already normalized. Since the BSDF is separable in θ and ϕ in this case, it is straightforward to determine the scatter angles from the random numbers:
The trajectory of the scattered ray s with respect to the surface normal n is given by:
where p and q represent the coordinate axes on the plane perpendicular to the surface normal. The trajectory of the scattered ray in (x,y,z) space is easily determined from the known components of (n,p,q) in this coordinate system:
An example: the Gaussian_XY scattering model
One of the example scattering DLLs provided with the OpticStudio installation is Gaussian_XY. The scattering distribution modeled by this DLL is a Gaussian about the projected specular ray. If p and q represent the coordinate axes on the scattering surface (i.e. perpendicular to the surface normal), then this distribution is represented analytically by:
As the probability is represented in Cartesian coordinates in this case, the corresponding integral is slightly simplified:
where Erf is the Error Function. Since P is separable, it is straightforward to assign values for the scatter direction based on random numbers:
Unfortunately the Error Function and its inverse are not built-in numeric functions in C. However, an alternate solution for this case exists:
This solution is based on the Box-Muller method for generating random numbers with a Gaussian distribution, as described in Chapter 7.2 of Numerical Recipes in C, Second Edition.1
The source code for this DLL (Gaussian_XY.c) is located in the Documents\Zemax\DLL\SurfaceScatter directory (assuming Opticstudio 17+). The main code starts with the Windows APIENTRY UserScatterDefinition. In the first part of this code, variables are defined, the random number generator is initialized, and the DLL reports to OpticStudio that a ray has been scattered and what the energy of that ray is (relative to the input energy available for scattering):
/* this DLL models a surface that has a Gaussian scattering
distribution about the projected specular ray */
int __declspec(dllexport) APIENTRY UserScatterDefinition(double *data)
{
double nx, ny, nz, sx, sy, sz, px, py, pz, qx, qy, qz, sig_x, sig_y;
double T, MAG, bo, random_number, random_val1, random_val2, bp, bq;
double tx, ty, tz;
const double PI = 4.0*atan(1.0);
/* we need some way of randomizing the random number generator */
srand((unsigned int) data[16]);
/* return transmission */
T = data[51];
if (T < 0.0 || T > 1.0) T = 1.0;
data[12] = T;
/* return flag to indicate that we scattered */
data[10] = 1.0;
Next, the coordinate axes p and q are determined such that they are perpendicular to the surface normal n, and such that the projected specular ray lies in the p-n plane:
/* find vectors p,q perpendicular to n, such that the specular ray lies along the p-n plane */
px = sx; // initialize p with values from the specular ray
py = sy;
pz = sz;
MAG = px*nx + py*ny + pz*nz;
if (MAG > 0.9999999)
{
/* specular is very close to the normal, any p will do */
if (fabs(nz) < 0.9)
{
px = 0.0;
py = 0.0;
pz = 1.0;
}
else
{
px = 1.0;
py = 0.0;
pz = 0.0;
}
}
/* this creates q normal to n and p */
CrossProduct(nx, ny, nz, px, py, pz, &qx, &qy, &qz);
Normalize(&qx, &qy, &qz);
/* this creates p normal to both q and n */
CrossProduct(nx, ny, nz, qx, qy, qz, &px, &py, &pz);
//Normalize(&px, &py, &pz); // not needed since n and q
// are orthonormal already
/* the specular ray, when projected onto the surface,lies along p */bo = px*sx + py*sy + pz*sz;
After the data for σp and σq are read in from the user input array, the scattered ray direction is determined using two random numbers:
get_another_scattered_ray:;
MAG = 2.0;
while (MAG > 1.0)
{
random_number = (double)rand()/(double)RAND_MAX;
if (random_number == 0.0) goto get_another_scattered_ray;
random_val1 = random_number;
random_number = (double)rand()/(double)RAND_MAX;
random_val2 = random_number;
bp = sig_x*sqrt(-1.0*log(random_val1))*cos(2.0*PI*random_val2); // Gaussian distribution along p
bq = sig_y*sqrt(-1.0*log(random_val1))*sin(2.0*PI*random_val2); // Gaussian distribution along q
tx = bo*px + bp*px + bq*qx; // t is the scattered ray
// vector projected onto
// the surface
ty = bo*py + bp*py + bq*qy;
tz = bo*pz + bp*pz + bq*qz;
MAG = tx*tx + ty*ty + tz*tz; // validate that scattered
// ray vector is within
// unit circle
}
MAG = sqrt(1.0 - MAG);
/* this is the scattered ray */
data[4] = tx + MAG*nx;
data[5] = ty + MAG*ny;
data[6] = tz + MAG*nz;
The random numbers are used to determine the scattering direction with respect to the projected specular ray according to the Box-Muller equations. The direction of the scattered vector in the p-q plane is then calculated using those values and the known direction of the projected specular ray in this plane. A check is performed to make sure that the scattered vector in the projected plane has a magnitude less than unity; if it does not, two new random numbers are drawn and the process is repeated. The actual trajectory of scattered ray in (x,y,z) space is finally determined by accounting for the component of scattering along the surface normal, which is calculated from the constraint that the magnitude of the scatter vector is unity. Once this calculation is complete, a check is performed to ensure that the ray did not scatter “beyond the horizon”, i.e. more than 90 degrees from the surface normal:
/*
NOW, make sure we didn't scatter in such a way that the ray is no more than 90 degrees from normal!
*/
if (nx*data[4] + ny*data[5] + nz*data[6] < 0.0)
{
goto get_another_scattered_ray;
}
return 0;
}
If this condition is violated, the process is repeated from the beginning by drawing two new random numbers; otherwise the scattering calculation is complete for the given input ray.
A simple example (Gaussian_XY.zmx) which illustrates the use of this DLL is provided in the .ZIP file located in the "Article Attachments" link at the top of this article. For σp = σq = 0.5, the radiant intensity distribution looks like:
References
1. Numerical Recipes in C, The Art of Scientific Computing, Second Edition, 1997.