Light scattered from surfaces can contribute significantly to the image quality and overall performance of an optical system. In many cases scattering is a desired effect, e.g. when a source with uniform angular distribution is created by passing a light beam through a Lambertian diffuser. Surface scattering distributions which are built into ZEMAX include the Lambertian, Gaussian, and ABg models. More information on these distributions may be found in Chapter 12 of the ZEMAX manual.

There may be cases where it is necessary to model objects whose surface scattering properties do not follow any of the built-in ZEMAX distributions. In such cases user-defined scattering models may be employed. This is done by constructing a Windows Dynamic Link Library (DLL), as described in the Knowledge Base article “How to Compile a User-Defined Surface”. 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 can construct the necessary DLL for you at a very competitive price.

There are a few examples of user-defined surface scattering models which are provided with your ZEMAX installation. Both the C source code and corresponding DLL file for each model are located in the directory \ZEMAX\Objects\DLL\SurfaceScatter\. The majority of each code is similar, containing “boilerplate” programming used to link between the DLL and ZEMAX. 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.

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(q,f) describes the probability of the ray scattering into a range of angles dq (centered about q) and df (centered about f):



where q is the polar angle with respect to the surface normal n and f is the azimuthal angle, as shown above. The necessary integral, performed in polar coordinates, is then:

I = ∫∫ P(q,f) sinq dq df

This indefinite integral (which is still a function of q and f) 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 q and f (i.e. the integration goes from 0 to p/2 in q and from 0 to 2p in f). The normalized integral (I_N = I/TIS) is then a function of q and f which – when taking into account the constant of integration that arises from the indefinite integration P(q,f) – is constrained to be between 0 and 1.

For the case when P(q,f) is separable in q and f – i.e. P(q,f) = F(q)*G(f) – I_N can be separated into two terms:

I_N(q,f) = I_q(q)*I_f(f)

Values for the scattered angles are determined by setting each of the two random numbers equal to I_q and I_f, respectively, and then inverting the equations to solve for q and f:

RN1 = I_q(q); q = I_q^-1(RN1)
RN2 = I_f(f); f = I_f^-1(RN2)

If P(q,f) is not separable, then an intermediate step is necessary, but the basic idea is the same. I_q(q,f) is first calculated by integrating over f:

I_q(q,f) = ∫ P(q,f) df

This indefinite integral is normalized to the definite integral (= I_q(q)) calculated by performing the same integration over the full limits of f (0 to 2p). I and TIS are then calculated by performing the indefinite and definite integration of I_q(q), respectively:

I = ∫ I_q(q) sinq dq
I_N = I/TIS

The value for the polar angle is determined by setting the first random number equal to I_N and inverting the equation to solve for q:

RN1 = I_N(q); q = I_N^-1(RN1)

The value for the azimuthal angle is then determined by setting the second random number equal to I_q(q,f) for the known value of q, and inverting to solve for f:

RN2 = I_q(q,f); f = I_q^-1(I_N^-1(RN1),RN2)

For this case when P(q,f) is not separable, the order of integration is in fact not important. We could have chosen to instead integrate P(q,f) over q first, and then over f. The technique for relating the random numbers RN1 and RN2 to q and f would have been identical:

RN1 = I_N(f); f = I_N^-1(RN1)
RN2 = I_f(q,f); q = I_f^-1(I_N^-1(RN1),RN2)

Integration may therefore be performed in whichever order is most convenient for the given probability function.

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 cosq factor:

I = ∫∫ BSDF(q,f) sinq cosq dq df

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 (A_p = A*cosq). Otherwise, the procedure for determining the scatter angles q and f 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/p. I is then given by:

I = ∫∫ BSDF(q,f) sinq cosq dq df 
  = (1/p)*[∫ sinq cosq dq]*[∫ df]
  = sin²q*(f/2p)

The total integrated scatter is calculated by setting the limits of integration from 0 to p/2 in q and from 0 to 2p in f:

TIS = [sin²q (0:p/2)]*[f/2p (0:2p)] = [1-0]*[1-0] = 1

Therefore I is already normalized. Since the BSDF is separable in q and f in this case, it is straightforward to determine the scatter angles from the random numbers:

RN1 = sin²q; q = sin^-1(√RN1)
RN2 = f/2p; f = 2p*RN2

The trajectory of the scattered ray s with respect to the surface normal n is given by:

s_n = cosq
s_p = sinq*cosf
s_q = sinq*sinf

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:

s_x = s_n*n_x + s_p*p_x + s_q*q_x
s_y = s_n*n_y + s_p*p_y + s_q*q_y
s_z = s_n*n_z + s_p*p_z + s_q*q_z

One of the example scattering DLLs provided with the ZEMAX 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:

P(p,q) = (4/(p*s_p*s_q))*exp[-((p/s_p)² + (q/s_q)²)]

As the probability is represented in Cartesian coordinates in this case, the corresponding integral is slightly simplified:

I = ∫∫ P(p,q) dp dq = Erf(p/s_p)*Erf(q/s_q)

where Erf is the Error Function. Since P is separable, it is straightforward to assign values for the scatter direction based on random numbers:

RN1 = Erf(p/s_p); p = s_p*Erf^-1(RN1)
RN2 = Erf(q/s_q); q = s_q*Erf^-1(RN2)

Unfortunately the Error Function and its inverse are not built-in numeric functions in C. However, an alternate solution for this case exists:

p = s_p*√(-ln(RV1))*cos(2p*RV2)
q = s_q*√(-ln(RV1))*sin(2p*RV2)

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¹.

The source code for this DLL (Gaussian_XY.c) is located in the \ZEMAX\Objects\DLL\SurfaceScatter\ directory. 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 ZEMAX 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 s_p and s_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 at the end of this article. For s_p = s_q = 0.5, the radiant intensity distribution looks like:

A number of built-in models exist in ZEMAX for describing scattering from a surface. However, there may be occasions when it is necessary to specify a scattering distribution not already provided. Scattering models can be easily added to ZEMAX by employing user-defined scattering functions, created as Windows DLL files. Information contained in the probability distribution or BSDF is used to perform scattering through a Monte Carlo technique. This technique simply requires that the distribution is integrable. If you require a user-defined scattering distribution, but do not wish to write the code yourself, we can construct the necessary DLL for you at a very competitive price.

References

¹Numerical Recipes in C, The Art of Scientific Computing, Second Edition, 1997.