Estimation of Lens Stray Light with Regard to the Incapacitation of Imaging Sensors

We present our efforts on estimating light scattering characteristics from commercial off-the-shelf (COTS) camera lenses in order to deduce thereof a set of generic scattering parameters valid for a specific lens class (double Gauss lenses). In previous investigations, we developed a simplified theoretical light scattering model to estimate the irradiance distribution in the focal plane of a camera lens. This theoretical model is based on a 3-parameter bidirectional scattering distribution function (BSDF), which describes light scattering from rough surfaces of the optical elements. Ordinarily, the three scatter parameters of the BSDF are not known for COTS camera lenses, which makes it necessary to assess them by own experiments. Besides the experimental setup and the measurement process, we present in detail the subsequent data exploitation. From measurements on seven COTS camera lenses, we deduced a generic set of scatter parameters. For a deeper analysis, the results of our measurements have also been compared with the output of an optical engineering software. Together with our theoretical model, now stray light calculations can be accomplished even then, when specific scatter parameters are not available from elsewhere. In addition, the light scattering analyses also allow considering the glare vulnerability of optical systems in terms of laser safety.


Introduction
Right after laser emission had been demonstrated the first time, the specific hazards of this new kind of light source became obvious [1,2]. Besides electrical risks associated with high voltage driven components, the main risk to individuals normally occurs through direct exposure of the eye to the laser beam. The worldwide efforts by a vast number of researchers to establish rules for the safe use of lasers led to the well-known laser safety standards like IEC 60825-1 or ANSI Z136.1 [3,4]. These standards provide quantities like maximum permissible exposure (MPE) limits, representing the highest level of irradiance or radiant exposure in order to enable a safe view into the laser beam for the human eye. The MPE depends on various parameters, like laser output power, wavelength, exposure time to the eye or pulse width and pulse repetition rate. Based on the MPE value and the beam divergence, the minimum distance between eye and laser source can be calculated, the so-called nominal ocular hazard distance (NOHD), below which the direct view into the laser beam is not safe. Recently, Williamson and McLin transferred the damage related MPE/NOHD concept to laser eye dazzle [5]. Equivalently to MPE and NOHD, they established the terms maximum dazzle exposure (MDE) and nominal ocular dazzle distance (NODD). Due to the increasing proliferation of hand-held high-power laser pointers and their often reported misuse, their work is of high importance regarding the evaluation of performance limitations in the execution of human tasks in cases of laser dazzle.
For imaging sensors, we encounter similar considerations with respect to sensor performance when they get dazzled. At our institute, we have been working for many years on the protection of imaging sensors from laser threats, comprising sensor hardening against laser damage and laser dazzle [6]. In a recent publication, one of the authors (G. R.) presented an approach for laser safety calculations for imaging sensors [7]. The approach transforms the above-mentioned quantities of laser safety calculations for the human eye to imaging sensors. The equivalent quantities derived for imaging sensors were called maximum permissible exposure for a sensor (MPE S ), nominal sensor hazard distance (NSeHD), maximum dazzle exposure for a sensor (MDE S ) and nominal sensor dazzle distance (NSeDD). We refer the reader to the publication of Ritt for details on the derivation of these quantities [7].
Very briefly, the derivation is based on the estimation of the radial irradiance distribution of laser light at (or near) the focal plane of a camera lens. This estimation considers diffraction of light at the lens aperture and scattering of light from the surfaces of the optical elements. The estimated irradiance distribution is then compared to threshold values for laser dazzle or laser damage of the sensor to estimate the already mentioned laser safety quantities. The primary goal of this approach was to establish closed-form expressions containing only basic operations and functions in order to calculate such quantities. Furthermore, to perform such calculations mainly the standard parameters of such devices are required, specified by the manufacturers of laser sources, camera lenses and cameras/imaging sensors. However, in general manufacturers do not provide information on the scattering characteristics of their lenses. That makes it difficult to estimate the fraction of light scattering in the irradiance distribution at the focal plane. Typically, three parameters are sufficient to describe light scattering of optical elements. In reference [7], results of initial scatter measurements were presented for commercial off-the-shelf (COTS) camera lenses. However, the main issue behind the problems touches the question of whether there is a strong variation of these parameters from camera lens to camera lens or may a single set of scatter parameters be sufficient to generally describe the behavior of light scattering inside COTS camera lenses. The answer has a significant impact on the ease applicability of our concept regarding laser safety calculations for imaging sensors.
Dedicated literature that would help to solve our tasks could not be found. There is a vast amount of publications on stray light. Typical values of scatter parameters of optical elements (mirrors or lenses) are stated in text books (e.g., reference [8]) or overview articles (e.g., reference [9]). Furthermore, there are many journal articles dedicated to stray light analysis of optical systems using optical engineering software like TracePro, FRED or ASAP. A larger number of them is related to astronomical optical systems, like mirror telescopes. There seems to be only a lower number of publications that is related to pure refractive optical systems, e.g., the camera lens of cellphones [10] or infrared imaging systems [11].
Usually, the stray light analysis using an optical engineering software requires the statement of scatter parameters or, alternatively, the default values implemented in the software can be taken. In all the publications we found, only standard values for the scatter parameters were applied. We did not find a single publication, where scatter parameters for the optical system to be simulated were measured before. Therefore, this situation forced us to make our own efforts to estimate scatter parameters of COTS camera lenses, as input to our theoretical model.
In this publication, we report on our measurements of the irradiance distribution in the focal plane of COTS camera lenses when illuminated with laser light. The results have been used to estimate the scatter parameters required by the theoretical model of reference [7]. In Section 2, a brief review on this theoretical model is given and explained by means of calculations of the focal plane irradiance of a standard camera. Sections 3 and 4 describe our experimental setup for the measurements and the camera lenses under test, respectively. Section 5 explains in detail the data analysis procedure for the estimation of the scatter parameters; the results are presented in Section 6. Finally, Section 7 treats with the simulation of stray light in camera lenses using the optical engineering software FRED in order to compare the theoretical model with the outcome of the FRED simulations.

Estimation of the Focal Plane Irradiance
The theoretical model of publication [7] assumes a scenario as depicted in Figure 1. A laser emits a beam with a Gaussian beam profile and illuminates a sensor consisting of a camera lens and an imaging sensor. In Figure 1, the camera lens is represented by a single lens, but is treated as an optical system consisting of several optical elements. In order to describe the irradiance distribution at the focal plane of the camera lens and the response of the imaging sensor to the laser radiation, we used the parameters listed in Table 1.

Estimation of the Focal Plane Irradiance
The theoretical model of publication [7] assumes a scenario as depicted in Figure 1. A laser emits a beam with a Gaussian beam profile and illuminates a sensor consisting of a camera lens and an imaging sensor. In Figure 1, the camera lens is represented by a single lens, but is treated as an optical system consisting of several optical elements. In order to describe the irradiance distribution at the focal plane of the camera lens and the response of the imaging sensor to the laser radiation, we used the parameters listed in Table 1.  (1 − exp (−2/ )) W Laser power entering the camera lens = m Laser spot size in the focal plane Spot size constant; see reference [12]   Truncation factor P in = P laser 1 − exp −2/ν 2 W Laser power entering the camera lens d spot = kλF m Laser spot size in the focal plane k Spot size constant; see reference [12] Sensors 2020, 20, 6308 4 of 42 Throughout this publication, we used the convention that the diameter of the laser beam is always the one at the position of the entrance aperture of the camera lens. Table 1 lists two definitions for the beam diameter. The first one, d 63 (m), is related to those points, where the irradiance has dropped to 1/e of the maximum irradiance. For Gaussian beams, 63 percent of the laser power is encircled within this diameter. When performing laser safety calculations, the use of the d 63 -diameter is mandatory. The second (more common) definition, d 86 (m), is related to those points, where the irradiance has dropped to 1/e 2 of the maximum irradiance. For Gaussian beams, 86 percent of the laser power is encircled within this diameter. The ratio of beam diameter d 86 to the diameter of the camera lens' aperture (more precisely it is the entrance pupil) d ap (m) is called the truncation factor ν and has a determining influence on the distribution of the laser light in the focal plane of the camera lens.
According to our theoretical model, the irradiance E fp (W/m 2 ) at the focal plane of a camera lens can be described by the sum of diffracted and scattered fractions of the incident light: This equation is identical to Equation (40) of reference [7] except that we here used the radial coordinate r (m) instead of the viewing angle Θ = r/ f (rad). On the right-hand side of Equation (1), the three quantities η d , E d (r) and E s (r) are given by Equations (39), (23) and (35a) of reference [7], respectively. There you also find details on the derivation of these equations. Briefly, η d describes the fraction of diffracted laser power. E d (r) is the irradiance distribution at the focal plane due to diffraction of a (truncated) Gaussian beam and E s (r) is the irradiance distribution at the focal plane due to light scattering from the surfaces of the various optical elements.
The term E d (r) is given by where E GA (r) represents the central lobe of the diffraction pattern approximated by a Gaussian curve. Since for typical camera/lens configurations, the aperture caused diffraction ring pattern is usually not fully resolved, E mean (r) describes the local mean irradiance of the wing parts of the diffraction pattern: The radial coordinate r pi , which separates the central lobe from the wing parts of the diffraction pattern, cannot be stated in the analytical form but has to be calculated numerically; see Equation (22) of reference [7].
How much of the incident radiation will be diffracted is determined by the term η d : where N ss is the number of scattering surfaces of the camera lens. TIS is the amount of total integrated scatter generated by one scattering surface [9]: Sensors 2020, 20, 6308

of 42
This quantity is a function of the three scatter parameters s, b (sr −1 ) and l (rad). These scatter parameters originate from the 3-parameter Harvey scatter model (see below) used to describe light scattering caused by the surface roughness of optical elements. The Harvey scatter model served Peterson to derive analytical equations to describe the distribution of scattered light in the focal plane of an optical system [13]. According to Peterson's model, the stray light irradiance at the focal plane E s (r) can be estimated by adding up the contributions E s,j (r) of the single scattering elements: where the contribution of a single scattering element is given by Here, NA is the numerical aperture, a ent (m) the radius of the beam at the first scattering element and a j (m) the radius of the beam at the jth scattering element. E ent (W/m 2 ) is the entering irradiance and b 0 is given by We applied the following simplification and approximation on Peterson's equations: Thus, we got rid of the dependence on the beam radius a j in Equation (8) and, after some transformations (for details see reference [7]), we obtained a simpler equation for the scattering caused irradiance distribution E s (r) in the focal plane of a camera lens due to scattering of laser light: where ν * is defined by Remark: Equation (11) corresponds to the form of the bidirectional scattering distribution function (BSDF) of the 3-parameter Harvey scatter model [8], expressed by when applying an incidence angle of Θ 0 = 0 and using the small angle approximation sin(Θ) ≈ Θ.
Choosing three different sets of scatter parameters b 0 , s and l, Figure 2 shows their influence on the BSDF. At low scatter angles Θ, the scattering signal stays constant and decreases with increasing values of Θ. In a double-logarithmic plot, the decay appears linear. The meaning of the three scatter parameters is depicted using the yellow curve in in Figure 2. Scatter parameter l indicates the scatter angle where the BSDF changes from the constant region to the decreasing part of the curve with slope s. In a double-logarithmic plot and for scatter angles Θ l, the scatter parameter s represents the slope of the BSDF. Scatter parameters b is the value of the BSDF for sin(Θ) − sin(Θ 0 ) = 0.01 and scatter parameter b 0 describes the maximum BSDF value. Generally, the scatter parameters are material parameters with fixed values but they are supposed to change with wavelength. The wavelength scaling laws are given by [14] ( ) = ( ) with as the reference wavelength, for which the scatter parameters are known, and is the wavelength, for which the scatter parameters shall be calculated.
Based on the equations described here, we can see that our theoretical model to estimate the focal plane irradiance distribution for camera lenses depends mainly on standard parameters of the devices, besides the three scatter parameters , , and the truncation factor .

Estimation of the Camera Response
The output of a digital camera is an image consisting of a number of pixels comprising digital gray values. Since we want to compare such digital images with the theoretical model described before in the further course of this publication, we need to transfer these digital values to irradiance values (or the other way round). For a camera that has a linear response to the number of photons arriving during exposure time, the digital signal can be calculated according to the EMVA 1288 standard by [15] where (W/m 2 ) is the irradiance at the pixel, (m 2 ) the pixel area, (s) the camera's exposure time, ℎ = 6.626 ⋅ 10 Js the Planck constant, = 2.99792458 ⋅ 10 m/s the vacuum speed of light, (m) the wavelength, the quantum efficiency, (DN/e-) the overall system gain and . (DN) the dark signal.
As an example, Figure 3 shows a set of four curves for the camera signal and irradiance as a function of radial coordinate for different values of the truncation factor . The curves were calculated using Equations (1) and (15) with the following set of parameters:  Generally, the scatter parameters are material parameters with fixed values but they are supposed to change with wavelength. The wavelength scaling laws are given by [14] b with λ 0 as the reference wavelength, for which the scatter parameters are known, and λ is the wavelength, for which the scatter parameters shall be calculated.
Based on the equations described here, we can see that our theoretical model to estimate the focal plane irradiance distribution for camera lenses depends mainly on standard parameters of the devices, besides the three scatter parameters s, b, l and the truncation factor ν.

Estimation of the Camera Response
The output of a digital camera is an image consisting of a number of pixels comprising digital gray values. Since we want to compare such digital images with the theoretical model described before in the further course of this publication, we need to transfer these digital values to irradiance values (or the other way round). For a camera that has a linear response to the number of photons arriving during exposure time, the digital signal can be calculated according to the EMVA 1288 standard by [15] where E (W/m 2 ) is the irradiance at the pixel, A (m 2 ) the pixel area, t exp (s) the camera's exposure time, h = 6.626 · 10 −34 Js the Planck constant, c = 2.99792458 · 10 8 m/s the vacuum speed of light, λ (m) the wavelength, η the quantum efficiency, K (DN/e-) the overall system gain and µ y.dark (DN) the dark signal.
As an example, Figure 3 shows a set of four curves for the camera signal µ y and irradiance E as a function of radial coordinate r for different values of the truncation factor ν. The curves were calculated using Equations (1) and (15) with the following set of parameters: • Laser: P laser = 1 µW, λ = 532 nm, d 86 = 25 mm; • Camera: p = 5 µm, η = 0.6, t exp = 1 ms, K = 0.4 DN/e − , bd = 12 bit, µ y.dark = 0; presented in Section 6. The scatter parameter = 2 mrad corresponds to a radial coordinate of = • = 200 µm or / = 40 pixel, which is marked by a dashed vertical line in Figure 3. The set of different f-numbers results in a set of different truncation factors = /( / ) of 0.5, 1.0, 2.0 and 4.0. Since the emitted laser power was kept constant for the calculation, the incident power entering the camera is different for each plotted curve and was calculated using the equation given in Table 1. The camera signal as plotted in the graph does not correspond to a real camera. A real camera is only able to generate signals within a certain dynamic range, marked in the graph by a gray background covering radial coordinates ≥ 1 pixel and signals of ∈ [1 DN, 2 − 1 DN]. For our example camera, signals above the upper limit of 2 − 1 DN = 4095 DN would be restricted to this value, i.e., the camera would be saturated. Signals

Experimental Setup and Measurement Procedure
A scheme of our experimental setup to measure the irradiance at the focal plane of a camera lens is shown in Figure 4. As a light source, we used a Toptica iChrome MLE multi-wavelength laser source. This laser source offers four different laser wavelengths of 488 nm, 515 nm, 561 nm and 640 nm, which all were coupled into a common single-mode fiber. Depending on the wavelength, the output power at the fiber exit port ranged from 40 to 100 mW. The laser light was collimated using a reflective fiber collimator FC (Thorlabs RC08APC-P01). Attenuator A1 (neutral density filter Thorlabs NE40B-A) was used to set the maximum laser power to a value in the order of 1 µW. Subsequently, the light path was divided using beam splitter BS (Thorlabs CM1-BP1), sending the reflected part to the reference photodiode PD (Ophir PD300R-UV sensor head with power meter Ophir Vega). The transmitted light passed a second attenuator A2, consisting of a set of neutral density filters (Thorlabs NExxB, where xx relates to the optical density). In the further course, the laser beam passed a folding mirror FM and was expanded by a Keplerian telescope with magnification = 6.7 built by a focusing lens L (Thorlabs LA1484-A, = 300 mm) and an off-axis parabolic mirror OPM (Optical The wavelength-dependent parameters of the optical elements (splitting ratio of the beam splitter BS, transmittance of the attenuator A2 and lens L, reflectivity of folding mirror FM and off- For the values of the scatter parameters, we anticipated mostly the results of our measurements presented in Section 6. The scatter parameter l = 2 mrad corresponds to a radial coordinate of r l = l· f = 200 µm or r l /p = 40 pixel, which is marked by a dashed vertical line in Figure 3. The set of different f-numbers F results in a set of different truncation factors ν = d 86 /( f /F) of 0.5, 1.0, 2.0 and 4.0. Since the emitted laser power was kept constant for the calculation, the incident power P in entering the camera is different for each plotted curve and was calculated using the equation given in Table 1.
The camera signal µ y as plotted in the graph does not correspond to a real camera. A real camera is only able to generate signals within a certain dynamic range, marked in the graph by a gray background covering radial coordinates r px ≥ 1 pixel and signals of µ y ∈ 1 DN, 2 bd − 1 DN . For our example camera, signals µ y above the upper limit of 2 12 − 1 DN = 4095 DN would be restricted to this value, i.e., the camera would be saturated. Signals µ y < 1 DN could only be measured as an average signal of a multitude of pixels. For a real camera, the upper and lower usable limits are given by the saturation gray value µ y.sat (DN) and the absolute sensitivity threshold µ y.min (DN), respectively. These values are slightly different to the theoretical values of 4095 DN and 1 DN.

Experimental Setup and Measurement Procedure
A scheme of our experimental setup to measure the irradiance at the focal plane of a camera lens is shown in Figure 4. As a light source, we used a Toptica iChrome MLE multi-wavelength laser source. This laser source offers four different laser wavelengths of 488 nm, 515 nm, 561 nm and 640 nm, which all were coupled into a common single-mode fiber. Depending on the wavelength, the output power at the fiber exit port ranged from 40 to 100 mW. The laser light was collimated using a reflective fiber collimator FC (Thorlabs RC08APC-P01). Attenuator A1 (neutral density filter Thorlabs NE40B-A) was used to set the maximum laser power to a value in the order of 1 µW. Subsequently, the light path was divided using beam splitter BS (Thorlabs CM1-BP1), sending the reflected part to the reference photodiode PD (Ophir PD300R-UV sensor head with power meter Ophir Vega). The transmitted light passed a second attenuator A2, consisting of a set of neutral density filters (Thorlabs NExxB, where xx relates to the optical density). In the further course, the laser beam passed a folding mirror FM and was expanded by a Keplerian telescope with magnification M = 6.7 built by a focusing lens L (Thorlabs LA1484-A, f 1 = 300 mm) and an off-axis parabolic mirror OPM (Optical Surfaces Ltd. 037-0220, f 2 = 2000 mm). Finally, the collimated laser beam was sent to the camera lens CL under test. The laser beam diameters at the entrance of the camera lens were 21. 5 axis parabolic mirror OPM) were calibrated before the measurements. Therefore, we were able to calculate the power within the laser beam at the position of the camera lens CL using the reading of the reference photodiode PD. In order to measure the stray light irradiance at the focal plane of the camera lens, we used camera C (Allied Vision Mako G-419B NIR, Stadtroda, Germany) as a detector. The parameters of the camera are given in Table 2, which lists parameters as specified by the manufacturer [16] and parameters measured by ourselves. The latter were measured according to the EMVA 1288 standard; see Appendix B. From Figure 3 we can see that the central lobe of the diffraction pattern is within a radius of about 1-2 pixels of our camera. Therefore, we could not expect to resolve the central lobe with our camera-based experimental setup.

DN
Each camera lens under test (with attached camera C) was centered with respect to the optical axis defined by the laser beam. Thus, the center of the laser spot coincided mostly with the center of the detector. Since the tested camera lenses were different in length, each camera lens was positioned so that the front facet of the camera lens coincided with that position on the optical axis where the The wavelength-dependent parameters of the optical elements (splitting ratio of the beam splitter BS, transmittance of the attenuator A2 and lens L, reflectivity of folding mirror FM and off-axis parabolic mirror OPM) were calibrated before the measurements. Therefore, we were able to calculate the power within the laser beam at the position of the camera lens CL using the reading of the reference photodiode PD.
In order to measure the stray light irradiance at the focal plane of the camera lens, we used camera C (Allied Vision Mako G-419B NIR, Stadtroda, Germany) as a detector. The parameters of the camera are given in Table 2, which lists parameters as specified by the manufacturer [16] and parameters measured by ourselves. The latter were measured according to the EMVA 1288 standard; see Appendix B. From Figure 3 we can see that the central lobe of the diffraction pattern is within a radius of about 1-2 pixels of our camera. Therefore, we could not expect to resolve the central lobe with our camera-based experimental setup. Each camera lens under test (with attached camera C) was centered with respect to the optical axis defined by the laser beam. Thus, the center of the laser spot coincided mostly with the center of the detector. Since the tested camera lenses were different in length, each camera lens was positioned Sensors 2020, 20, 6308 9 of 42 so that the front facet of the camera lens coincided with that position on the optical axis where the laser beam diameter was measured. The camera and the laser were switched on 30-60 min before performing the measurements to ensure thermal equilibrium conditions. Figure 5 shows some photographs of the experimental setup. The complete experimental setup was covered by a housing to prevent from ambient light. We also made sure that the residual light in the laboratory (e.g., due to the computer monitor, emergency exit lights, etc.) had no noticeable effect on the camera signal. All experimental parameters were controlled using a computer, whereby the lens' f-number had to be set manually.
Sensors 2020, 20, x FOR PEER REVIEW 9 of 41 laser beam diameter was measured. The camera and the laser were switched on 30-60 min before performing the measurements to ensure thermal equilibrium conditions. Figure 5 shows some photographs of the experimental setup. The complete experimental setup was covered by a housing to prevent from ambient light. We also made sure that the residual light in the laboratory (e.g., due to the computer monitor, emergency exit lights, etc.) had no noticeable effect on the camera signal. All experimental parameters were controlled using a computer, whereby the lens' f-number had to be set manually. The main challenge in measuring the irradiance distribution of a focused laser beam is the high dynamic range of irradiance values, which have to be covered. At the center of the laser spot, the irradiance is quite high (10 4 W/m 2 in Figure 3), whereas the off-center stray light irradiance is quite low (10 −4 W/m 2 in Figure 3). Typically, the dynamic range of camera sensors is in the order of 60 dB, which means that the ratio of highest to lowest measureable irradiance corresponds to a factor of only 1000. This is by far not high enough to measure the irradiance within the whole area of the imaging sensor by capturing just one single image. Therefore, to gain an image with full intensity distribution, i.e., full dynamic range, we had to acquire a number of camera images based on different combinations of laser power and camera exposure time. For the further course of the publication, we introduced a naming convention for the different measurement steps, which are explained below in more detail; an overview is given in Table 3.
For each camera lens under test, we started the experiments by setting a specific f-number . Using a selected setting of camera exposure time and laser power , in the first step, we performed the image acquisition process. Dependent on the exposure time and laser power used, some parts of such an image may be overexposed or too noisy. Thus, such a single image will only deliver a part of the complete radial irradiance profile with a linear signal response. In order to get the desired irradiance information for almost the complete area of the imaging sensor, this image acquisition process was repeated seven times for different combinations of exposure time and laser power (see Table 4) to obtain a total of eight image acquisitions. We call this process of eight image acquisitions a measurement. Then, such a measurement performed for each of the four available laser wavelengths, results in a measurement series. By repeating such a measurement series for each labeled f-number of a camera lens, we obtained a complete data set for this specific camera lens. The output of this procedure accomplished for all the mentioned camera lenses (see Section 4) eventually formed a data ensemble.
Besides the image acquisitions with laser illumination, we also acquired images without illumination (dark frames), necessary for the later data analysis. The dark frames were recorded according to the exposure times listed in Table 4. The main challenge in measuring the irradiance distribution of a focused laser beam is the high dynamic range of irradiance values, which have to be covered. At the center of the laser spot, the irradiance is quite high (10 4 W/m 2 in Figure 3), whereas the off-center stray light irradiance is quite low (10 −4 W/m 2 in Figure 3). Typically, the dynamic range of camera sensors is in the order of 60 dB, which means that the ratio of highest to lowest measureable irradiance corresponds to a factor of only 1000. This is by far not high enough to measure the irradiance within the whole area of the imaging sensor by capturing just one single image. Therefore, to gain an image with full intensity distribution, i.e., full dynamic range, we had to acquire a number of camera images based on different combinations of laser power and camera exposure time. For the further course of the publication, we introduced a naming convention for the different measurement steps, which are explained below in more detail; an overview is given in Table 3.
For each camera lens under test, we started the experiments by setting a specific f-number F. Using a selected setting of camera exposure time t exp and laser power P laser , in the first step, we performed the image acquisition process. Dependent on the exposure time and laser power used, some parts of such an image may be overexposed or too noisy. Thus, such a single image will only deliver a part of the complete radial irradiance profile with a linear signal response. In order to get the desired irradiance information for almost the complete area of the imaging sensor, this image acquisition process was repeated seven times for different combinations of exposure time t exp and laser power P laser (see Table 4) to obtain a total of eight image acquisitions. We call this process of eight image acquisitions a measurement. Then, such a measurement performed for each of the four available laser wavelengths, results in a measurement series. By repeating such a measurement series for each labeled f-number F of a camera lens, we obtained a complete data set for this specific camera lens. The output of this procedure accomplished for all the mentioned camera lenses (see Section 4) eventually formed a data ensemble. Besides the image acquisitions with laser illumination, we also acquired images without illumination (dark frames), necessary for the later data analysis. The dark frames were recorded according to the exposure times listed in Table 4. Figure 6 shows some example images for the camera lens Edmund Optics 86410, acquired at the laser wavelength of λ = 488 nm and for an f-number of F = 2.8. The red labels correspond to the setting number indicated in Table 4. Please note, that the individual images do not show the complete camera image comprising 2048 pixels × 2048 pixels, but only the central section. The individual image sections also differed in size depending on the exposed area; see the white annotations within the images.    Figure 6 shows some example images for the camera lens Edmund Optics 86410, acquired at the laser wavelength of = 488 nm and for an f-number of = 2.8. The red labels correspond to the setting number indicated in Table 4. Please note, that the individual images do not show the complete camera image comprising 2048 pixels × 2048 pixels, but only the central section. The individual image sections also differed in size depending on the exposed area; see the white annotations within the images.  Table 4. Camera lens: Edmund Optics 86410, experimental parameters: λ = 488 nm, F = 2.8. Table 5 lists the lenses used for the experiments and their specifications. Intentionally, low-priced and higher-priced camera lenses were chosen: seven commercial off-the-shelf (COTS) camera lenses and one COTS achromatic doublet lens. The lenses differed in their values of focal length, f-number and number of optical elements. The set of f-numbers used for the measurements varied for each camera lens. We only selected those settings that were labeled on the aperture ring. One camera lens (Edmund Optics 54690) had no labels at the aperture ring. Here, we only used the minimum and the maximum f-number (largest und smallest aperture) for the experiments. In case of the achromatic doublet Thorlabs AC254-050-A, we used an external iris (Thorlabs SM2D25D) directly in front of the lens to obtain different settings for the f-number. Figure 7 shows photographs of all lenses and in addition, details of their aperture rings. For the achromatic doublet, the image shows the complete assembly with tube housing and an external iris. The labels correspond to the setting numbers of Table 4. Camera lens: Edmund Optics 86410, experimental parameters: = 488 nm, = 2.8. Table 5 lists the lenses used for the experiments and their specifications. Intentionally, lowpriced and higher-priced camera lenses were chosen: seven commercial off-the-shelf (COTS) camera lenses and one COTS achromatic doublet lens. The lenses differed in their values of focal length, fnumber and number of optical elements.

Camera Lenses Tested
The set of f-numbers used for the measurements varied for each camera lens. We only selected those settings that were labeled on the aperture ring. One camera lens (Edmund Optics 54690) had no labels at the aperture ring. Here, we only used the minimum and the maximum f-number (largest und smallest aperture) for the experiments. In case of the achromatic doublet Thorlabs AC254-050-A, we used an external iris (Thorlabs SM2D25D) directly in front of the lens to obtain different settings for the f-number. Figure 7 shows photographs of all lenses and in addition, details of their aperture rings. For the achromatic doublet, the image shows the complete assembly with tube housing and an external iris.    Table 6 lists the settings of the f-numbers used for the experiments and the corresponding mean truncation factors. Since the laser beam diameter slightly varied with the wavelength, the truncation factor also varied slightly with the wavelength and therefore only the mean truncation factors v are given. Table 6. f-number settings used for the measurement series and their corresponding mean truncation factor.

Data Analysis
The aim of our investigations was to estimate for each camera lens the scatter parameters s, b and l, in order to generate a generic set of scatter parameters. In the course of our experiments, we analyzed eight different lenses (seven camera lenses and one achromatic doublet, see Section 4) applying four different laser wavelengths at various integration times and laser powers, while the laser beam diameter was kept constant throughout. For each lens, the data set contained 1-7 measurement series, depending on the number of labeled settings of the lens' aperture ring. The data ensemble of eight data sets contained in total 43 measurement series (see Table 6), i.e., 43 × 4 = 172 measurements or 43 × 4 × 8 = 1376 camera images, to be analyzed. Since the amount of data was quite huge, we used automated analysis software to derive the scatter parameters. Briefly, this included the following steps: 1.
Irradiance profile generation: For each acquired image, an irradiance profile was assessed from the image data.

2.
Profile stitching: Due to the limited dynamic range of the camera, the merging of the eight radial irradiance profiles of a measurement was necessary to get a complete radial irradiance profile for the specific experimental set-up (camera lens under test, f-number F and laser wavelength λ).

3.
Curve fitting: Fitting of a theoretical curve to the radial irradiance profiles of a measurement series (comprising the measurements with all four laser wavelengths) using the scatter parameters s, b and l as fitting parameters.
Now, we will describe in detail the above-mentioned analysis process. As an example, we used part of the data of the camera lens Edmund Optics 86410. Example images acquired with this camera lens for a laser wavelength of λ = 488 nm and an f-number F = 2.8 are shown in Figure 6. Subsequent to the data analysis process for all camera lenses, we performed a statistical analysis of the scatter parameters. These results are presented in Section 6.

Dark Frame Correction
Before analyzing the image data, we performed dark frame corrections on the laser-illuminated images, i.e., the dark frames were subtracted from the laser-illuminated images to remove the dark Sensors 2020, 20, 6308 13 of 42 signal. The residual signal of the camera pixels is then proportional to the incident irradiance; see the camera's linearity in Appendix B.

Estimation of the Center of the Laser Spot
As a prerequisite for the generation of the irradiance profile, the center of the laser spot within the camera images had to be determined first. Although the alignment of the camera (with attached camera lens) ensured that the position of the laser spot was always close to the center of the imaging sensor (see Section 3) the slight variations regarding the exact center could not be neglected.
To find the center of the laser spot the following procedure was applied: we chose the three images acquired at the lowest laser power, using settings no. 1-3 of Table 4. From each of these three images, we extracted the central part of 100 pixels × 100 pixels and calculated a mean image. Subsequently, the pixel with a maximum signal within this mean image (coordinates: column x c , row y c ) was identified and defined as the center of the laser spot.
For the camera lens Edmund Optics 86410 (F = 2.8) this procedure is depicted in Figure 8 for the two laser wavelengths of 488 nm and 515 nm. In each frame, a red cross marks the estimated center of the laser spot.

Estimation of the Center of the Laser Spot
As a prerequisite for the generation of the irradiance profile, the center of the laser spot within the camera images had to be determined first. Although the alignment of the camera (with attached camera lens) ensured that the position of the laser spot was always close to the center of the imaging sensor (see Section 3) the slight variations regarding the exact center could not be neglected.
To find the center of the laser spot the following procedure was applied: we chose the three images acquired at the lowest laser power, using settings no. 1-3 of Table 4. From each of these three images, we extracted the central part of 100 pixels × 100 pixels and calculated a mean image. Subsequently, the pixel with a maximum signal within this mean image (coordinates: column , row ) was identified and defined as the center of the laser spot.
For the camera lens Edmund Optics 86410 ( = 2.8) this procedure is depicted in Figure 8 for the two laser wavelengths of 488 nm and 515 nm. In each frame, a red cross marks the estimated center of the laser spot.  Table 4 and the corresponding mean image. A red cross marks the center of the laser spot (estimated from the mean image). The red crosses in the images #1-#3 are derived from the mean image.

Estimation of the Radial Irradiance Profile
In order to find the lens-generated radial irradiance profile in the focal plane, the pixel values for each occurring radial distance to the center of the laser sport were averaged. The principle of the process is illustrated in Figure 9. It shows a complete image of 2048 pixels × 2048 pixels, as taken with setting no. 8 in Figure 6. A red arrow indicates the radial coordinate axis. The orange arrow depicts the circular averaging process for a specific radial coordinate.  Table 4 and the corresponding mean image. A red cross marks the center of the laser spot (estimated from the mean image). The red crosses in the images #1-#3 are derived from the mean image.

Estimation of the Radial Irradiance Profile
In order to find the lens-generated radial irradiance profile in the focal plane, the pixel values for each occurring radial distance r to the center of the laser sport were averaged. The principle of the process is illustrated in Figure 9. It shows a complete image of 2048 pixels × 2048 pixels, as taken with setting no. 8 in Figure 6. A red arrow indicates the radial coordinate axis. The orange arrow depicts the circular averaging process for a specific radial coordinate.
The averaging process first demands to calculate the radial distance r px (pixel) of each pixel in the frame regarding the center of the laser spot (x c , y c ): Sensors 2020, 20, x FOR PEER REVIEW 14 of 41 Figure 9. Illustration of the averaging process. The signals of all those pixels having the same distance (red coordinate axis) to the laser spot center are averaged (depicted by the circular arrow in orange color) in order to find the radial irradiance distribution. Colored disks, overlaid to the camera image, depict areas where different rounding precisions were applied. The camera image was taken using the lens Edmund Optics 86410 ( = 488 nm, = 2.8, setting no. 8).
The averaging process first demands to calculate the radial distance (pixel) of each pixel in the frame regarding the center of the laser spot ( , ): As one can easily see, for each occurring value of distance there are only eight pixels with exactly the same distance to the center of the laser spot. This would have led to a set of data points for the irradiance profile with a large number of values for the independent variable , but a low statistical basis for the dependent variable . Therefore, before averaging took place the values of were rounded. The precision for rounding was adapted to different ranges of the radial coordinate:
The areas with different rounding precision are marked in Figure 9 by colored disks, overlaid to the camera image. The rationale behind this choice of these areas will become obvious in Section 5.3.1.
The rounding procedure is illustrated in Figure 10 for better understanding. The squares represent individual pixels out of the 2048 pixels × 2048 pixels of the imaging sensor, whereby the three pixel blocks represent three distant regions where the rounding precision was adapted to an individually chosen value, i.e., rounding precision of 0.5, 1.0 and 2.0. The coordinate axes indicate the / -coordinates of the pixels regarding the laser spot center ( , ), which is highlighted by the red frame in the left pixel block. Inside the squares, the black numbers show the radial coordinate , exactly calculated to six decimal places while the red numbers show the corresponding rounded value of , used for later analysis according to the above-described procedure. Figure 9. Illustration of the averaging process. The signals of all those pixels having the same distance r px (red coordinate axis) to the laser spot center are averaged (depicted by the circular arrow in orange color) in order to find the radial irradiance distribution. Colored disks, overlaid to the camera image, depict areas where different rounding precisions were applied. The camera image was taken using the lens Edmund Optics 86410 (λ = 488 nm, F = 2.8, setting no. 8).
As one can easily see, for each occurring value of distance r px there are only eight pixels with exactly the same distance r px to the center of the laser spot. This would have led to a set of data points for the irradiance profile with a large number of values for the independent variable r px , but a low statistical basis for the dependent variable µ y . Therefore, before averaging took place the values of r px were rounded. The precision for rounding was adapted to different ranges of the radial coordinate:
The areas with different rounding precision are marked in Figure 9 by colored disks, overlaid to the camera image. The rationale behind this choice of these areas will become obvious in Section 5.3.1.
The rounding procedure is illustrated in Figure 10 for better understanding. The squares represent individual pixels out of the 2048 pixels × 2048 pixels of the imaging sensor, whereby the three pixel blocks represent three distant regions where the rounding precision was adapted to an individually chosen value, i.e., rounding precision of 0.5, 1.0 and 2.0. The coordinate axes indicate the x/y-coordinates of the pixels regarding the laser spot center (x c , y c ), which is highlighted by the red frame in the left pixel block. Inside the squares, the black numbers show the radial coordinate r px , exactly calculated to six decimal places while the red numbers show the corresponding rounded value of r px , used for later analysis according to the above-described procedure.
For an acquired image, the successive circular averaging processes resulted in a radial irradiance profile. According to this procedure, the image data of Figure 6 was treated correspondingly and the respective radial irradiance profiles are shown in Figure 11a. The results shown there, we denote as raw data, since the curves contain also measurement values not usable for the further processing. For an acquired image, the successive circular averaging processes resulted in a radial irradiance profile. According to this procedure, the image data of Figure 6 was treated correspondingly and the respective radial irradiance profiles are shown in Figure 11a. The results shown there, we denote as raw data, since the curves contain also measurement values not usable for the further processing. From Figure 11a, we can clearly recognize the saturated parts of the various irradiance profiles (course of the curves is horizontal) and parts that look noisy (strong signal scattering). Therefore, we filtered the data and kept only those values that belong to the slope of the curves: 1. We dismissed those radial coordinates for which a. More than 10 percent of the pixels have gray values larger than the saturation gray value 2. For the signals of those radial coordinates that were not dismissed, we calculated the average using only those pixel values that were within the limits [ . , . ], i.e., we calculated a trimmed mean. For an acquired image, the successive circular averaging processes resulted in a radial irradiance profile. According to this procedure, the image data of Figure 6 was treated correspondingly and the respective radial irradiance profiles are shown in Figure 11a. The results shown there, we denote as raw data, since the curves contain also measurement values not usable for the further processing. From Figure 11a, we can clearly recognize the saturated parts of the various irradiance profiles (course of the curves is horizontal) and parts that look noisy (strong signal scattering). Therefore, we filtered the data and kept only those values that belong to the slope of the curves: 1. We dismissed those radial coordinates for which a. More than 10 percent of the pixels have gray values larger than the saturation gray value 2. For the signals of those radial coordinates that were not dismissed, we calculated the average using only those pixel values that were within the limits [ . , . ], i.e., we calculated a trimmed mean. From Figure 11a, we can clearly recognize the saturated parts of the various irradiance profiles (course of the curves is horizontal) and parts that look noisy (strong signal scattering). Therefore, we filtered the data and kept only those values that belong to the slope of the curves:

1.
We dismissed those radial coordinates for which a.
More than 10 percent of the pixels have gray values larger than the saturation gray value µ y.sat (DN) or; b. More than 10 percent of the pixels have gray values lower than the absolute sensitivity threshold µ y.min (DN).

2.
For the signals of those radial coordinates that were not dismissed, we calculated the average using only those pixel values that were within the limits µ y.min , µ y.sat , i.e., we calculated a trimmed mean.
The saturation gray value of our camera was estimated to be µ y.sat = 3861; see Appendix B. The absolute sensitivity threshold µ y.min (related to the digital gray value) was calculated from the corresponding value µ e.min (related to the number of photoelectrons) by where K (DN/e-) is the overall system gain. For our camera, the absolute sensitivity threshold was µ e.min = 14.1 e − and the overall system gain was K = 0.399 DN/e − (see Table 2). Thus, for the sensitivity threshold we get µ y.min = 0.399 DN/e − · 14.1 e − ≈ 6 DN. Furthermore, we only kept that part of the irradiance profile, where the radial coordinates do not exceed the edges of the camera image in any direction. This means the camera signals in the corners of the image were discarded. The maximum value of r px is depicted in Figure 9 by the outer edge of the blue disk for a rounding precision of 2.0. For a perfect alignment of the laser spot to the exact center of the imaging sensor, the maximum value of r px would correspond to 1023 pixels. The residual data after the complete filtering process is shown in Figure 11b.

Profile Stitching
To get the complete radial irradiance profile, the individual profile sections were stitched together. Since they were derived for different settings of the camera's exposure time and the optical density of attenuator A2, we had to correct the values accordingly. For this, we scaled the values by a factor where i is the setting number according to Table 4, P i the laser power during image acquisition for setting number i and P = 8 i=1 P i /8 is the mean power during the whole measurement. The individual radial irradiance profile sections of Figure 11b scaled in this way are shown in Figure 12. Since the individual profiles may overlap regarding the radial coordinate r, we additionally calculated the mean of those values having the same radial coordinate. The resulting complete irradiance profile is shown in Figure 12 by the black solid line and was used for the further analysis process. It would correspond to a profile a sensor would provide when there were no limitations in its dynamic range and using an exposure time of 1 µs. The saturation gray value of our camera was estimated to be . = 3861; see Appendix B. The absolute sensitivity threshold .
(related to the digital gray value) was calculated from the corresponding value .
(related to the number of photoelectrons) by where (DN/e-) is the overall system gain. For our camera, the absolute sensitivity threshold was . = 14.1 e and the overall system gain was = 0.399 DN/e (see Table 2). Thus, for the sensitivity threshold we get . = 0.399 DN/e ⋅ 14.1 e ≈ 6 DN.
Furthermore, we only kept that part of the irradiance profile, where the radial coordinates do not exceed the edges of the camera image in any direction. This means the camera signals in the corners of the image were discarded. The maximum value of is depicted in Figure 9 by the outer edge of the blue disk for a rounding precision of 2.0. For a perfect alignment of the laser spot to the exact center of the imaging sensor, the maximum value of would correspond to 1023 pixels. The residual data after the complete filtering process is shown in Figure 11b.

Profile Stitching
To get the complete radial irradiance profile, the individual profile sections were stitched together. Since they were derived for different settings of the camera's exposure time and the optical density of attenuator A2, we had to correct the values accordingly. For this, we scaled the values by a factor , 1 µ ⋅ 10 , ⋅ , where is the setting number according to Table 4, the laser power during image acquisition for setting number and = ∑ /8 is the mean power during the whole measurement. The individual radial irradiance profile sections of Figure 11b scaled in this way are shown in Figure 12.
Since the individual profiles may overlap regarding the radial coordinate , we additionally calculated the mean of those values having the same radial coordinate. The resulting complete irradiance profile is shown in Figure 12 by the black solid line and was used for the further analysis process. It would correspond to a profile a sensor would provide when there were no limitations in its dynamic range and using an exposure time of 1 µs.

Curve Fitting
Subsequent to the estimation of the complete radial irradiance profiles, different curves based on different theoretical models were fitted to the data in order to find the desired scatter parameters. Here, we described our approach and the rationales behind it.

Curve Fitting
Subsequent to the estimation of the complete radial irradiance profiles, different curves based on different theoretical models were fitted to the data in order to find the desired scatter parameters. Here, we described our approach and the rationales behind it.
For fitting the data, we used three different theoretical models, which we will discuss at a later stage in detail: All models have in common that they describe the spatial distribution of the irradiance E in the focal plane of a camera lens as function of the radial coordinate r (m) and assuming rotational symmetry. That makes it necessary to adapt the equations of the corresponding models to the irradiance profiles that represent the camera signal µ y in units of digital numbers (DN) as function of the radial coordinate r px (pixel). For this, we related the radial coordinate r in our theoretical models to the radial coordinate r px by where p (m) is the pixel size of the camera (p = 5.5 µm). Furthermore, the irradiance values had to be transformed into camera signals using Equation (15) with µ y.dark = 0, because of the dark frame correction: Due to the profile stitching process, see Section 5.2, we used a fixed exposure time of t exp = 1 µs.

Model M1: Our Original Model
In a first step, we used our original theoretical model M1 to simulate the focal plane irradiance as described by Equation (1): The scatter parameters s, b and l of the term E s (r) were used as fit parameters. The fitting ranges for these parameters were [−3.5 ≤ s ≤ −0.5], [0.01 ≤ b ≤ 100] and 5 · 10 −4 ≤ l ≤ 0.1 .
In principle, we could perform a fit on the data of each measurement (single wavelength). Alternatively, we could perform a fit on the complete data of the measurement series (comprising all four laser wavelengths) by including the wavelength-scaling laws for the scatter parameters given by Equation (14). We recognized that the latter method leads to a more robust fit and, thus, used this method for the data analyses. The outcome of the fitting process was always related to a reference wavelength of 550 nm.
As an example for the curve fitting according to model M1, Figure 13 shows the derived irradiance profiles for the camera lens Edmund Optic 86410 for two different values of the f-number; Figure 13a Figure 12 (λ = 488 nm) can be found as blue data points in Figure 13a. The black lines show the model curves after the fitting process. The value of the scatter parameter l is indicated by the vertical lines; the color coding of these lines corresponds to that of the data points.
By the example of Figure 13, we can see that the radial range of 10 pixel < r px < 200 pixel represents a transition range, which divides the region where mainly scatter dominates (r px > 200 pixel) from the region where mainly diffraction (and aberration) dominates (r px < 10 pixel). This transition is much more pronounced for smaller values of the truncation factor ν. For very high values of the Sensors 2020, 20, 6308 18 of 42 f-number, the transition may even not be clearly visible. These observations apply to all the camera lenses and are the reason for our choice to apply different values of rounding precisions for the profile generation process described in Section 5.1. Most of the data is located in the region r px > 200 pixel. By our choice of the rounding precisions, we were able to adjust the number of data points in such a way that the region r px > 200 pixel does not dominate the fitting result and that the characteristics of the curve for r px < 200 pixel also influence the outcome of the curve fit. By the example of Figure 13, we can see that the radial range of 10 pixel < < 200 pixel represents a transition range, which divides the region where mainly scatter dominates ( > 200 pixel) from the region where mainly diffraction (and aberration) dominates ( < 10 pixel). This transition is much more pronounced for smaller values of the truncation factor . For very high values of the f-number, the transition may even not be clearly visible. These observations apply to all the camera lenses and are the reason for our choice to apply different values of rounding precisions for the profile generation process described in Section 5.1. Most of the data is located in the region > 200 pixel. By our choice of the rounding precisions, we were able to adjust the number of data points in such a way that the region > 200 pixel does not dominate the fitting result and that the characteristics of the curve for < 200 pixel also influence the outcome of the curve fit.
Furthermore, we can see that there is a discrepancy between the irradiance profile and the model for very small radial coordinates ( < 10 pixel). We attributed this to the fact that model M1, described by Equation (1), did not include aberrations. However, aberrations will affect the center of the laser spot. The observed discrepancy is stronger for small values of truncation factor . This discrepancy leads to strong variations of the fitting results for different values of the f-number; see, e.g., the different location of the vertical lines in Figure 13a,b. For this reason, we adapted our analysis process and extended model M1 by an additional term (then named model M2) in order to minimize the deviations between the theoretical model and the measured data in the vicinity of the center of the laser spot.

Model M2: The Auxiliary Model
By extending model M1 by an empirical additional term, we got model M2, which may be seen as an approach to account for the previous discrepancies between model M1 and the measured data at the center of the laser spot: The additional term (red part in Equation (21)) extends the diffraction term ( ) for aberrations by multiplying it with a Gaussian function. Mathematically, the two parameters and describe the amplitude and the width of the Gaussian function, respectively, and were used as additional fitting parameters. The fitting parameters and were kept constant for all laser Furthermore, we can see that there is a discrepancy between the irradiance profile and the model for very small radial coordinates (r px < 10 pixel). We attributed this to the fact that model M1, described by Equation (1), did not include aberrations. However, aberrations will affect the center of the laser spot. The observed discrepancy is stronger for small values of truncation factor ν. This discrepancy leads to strong variations of the fitting results for different values of the f-number; see, e.g., the different location of the vertical lines in Figure 13a,b.
For this reason, we adapted our analysis process and extended model M1 by an additional term (then named model M2) in order to minimize the deviations between the theoretical model and the measured data in the vicinity of the center of the laser spot.

Model M2: The Auxiliary Model
By extending model M1 by an empirical additional term, we got model M2, which may be seen as an approach to account for the previous discrepancies between model M1 and the measured data at the center of the laser spot: The additional term (red part in Equation (21)) extends the diffraction term η d E d (r) for aberrations by multiplying it with a Gaussian function. Mathematically, the two parameters p 1 and p 2 describe the amplitude and the width of the Gaussian function, respectively, and were used as additional fitting parameters. The fitting parameters p 1 and p 2 were kept constant for all laser wavelengths. The auxiliary model M2 is able to describe the measurement results much better near the center of the laser spot as compared to the original model M1. Figure 14 presents the results achieved with model M2 for the same examples as given in Figure 13. Since we had no real physical explanation for the new parameters and in the auxiliary model, we just used the fitting parameter to set a lower limit of the radial coordinate to perform alternative curve fits with model M1, but now based on an adjusted radial coordinate range. We defined the range where the measurement values deviated from the theoretical model M1 as [0, 2 ] and set the minimum value of the radial coordinate to = 2 for the second fitting process with model M1. Results are presented in Figure 15 for the same examples as used before, where the excluded data is indicated in the plots by a gray background. We can see that now the location of the vertical lines is comparable for the two different f-numbers.

Model M3: The Simplified Model
The derivation of the scatter component ( ) of our original model M1 required taking the diameter of the incident light beam at the optics into account [7]. Two cases had to be distinguished, see Equation (10): (i) For -beam diameters smaller than the lens aperture we used as the beam diameter for the calculations. (ii) For beam diameters larger than the lens aperture , the beam diameter was set to the value of the lens aperture . This definition led to the factor * of Equation (12). Now, as an assumption, we set ⋆ ≔ 1. This means that the incident light beam is assumed to fill the complete lens aperture. For Gaussian laser beams, this assumption is not false since the wings of a Gaussian beam profile extends, at least in theory, to infinity. Such an assumption Since we had no real physical explanation for the new parameters p 1 and p 2 in the auxiliary model, we just used the fitting parameter p 2 to set a lower limit of the radial coordinate r to perform alternative curve fits with model M1, but now based on an adjusted radial coordinate range. We defined the range where the measurement values deviated from the theoretical model M1 as [0, 2p 2 ] and set the minimum value of the radial coordinate to r min = 2p 2 for the second fitting process with model M1. Results are presented in Figure 15 for the same examples as used before, where the excluded data is indicated in the plots by a gray background. We can see that now the location of the vertical lines is comparable for the two different f-numbers. Since we had no real physical explanation for the new parameters and in the auxiliary model, we just used the fitting parameter to set a lower limit of the radial coordinate to perform alternative curve fits with model M1, but now based on an adjusted radial coordinate range. We defined the range where the measurement values deviated from the theoretical model M1 as [0, 2 ] and set the minimum value of the radial coordinate to = 2 for the second fitting process with model M1. Results are presented in Figure 15 for the same examples as used before, where the excluded data is indicated in the plots by a gray background. We can see that now the location of the vertical lines is comparable for the two different f-numbers.

Model M3: The Simplified Model
The derivation of the scatter component ( ) of our original model M1 required taking the diameter of the incident light beam at the optics into account [7]. Two cases had to be distinguished, see Equation (10): (i) For -beam diameters smaller than the lens aperture we used as the beam diameter for the calculations. (ii) For beam diameters larger than the lens aperture , the beam diameter was set to the value of the lens aperture . This definition led to the factor * of Equation (12). Now, as an assumption, we set ⋆ ≔ 1. This means that the incident light beam is assumed to fill the complete lens aperture. For Gaussian laser beams, this assumption is not false since the wings of a Gaussian beam profile extends, at least in theory, to infinity. Such an assumption

Model M3: The Simplified Model
The derivation of the scatter component E s (r) of our original model M1 required taking the diameter of the incident light beam at the optics into account [7]. Two cases had to be distinguished, see Equation (10): (i) For d 63 -beam diameters smaller than the lens aperture d ap we used d 63 as the beam diameter for the calculations. (ii) For beam diameters d 63 larger than the lens aperture d ap , the beam diameter was set to the value of the lens aperture d ap . This definition led to the factor ν * of Equation (12). Now, as an assumption, we set ν := 1. This means that the incident light beam is assumed to fill the complete lens aperture. For Gaussian laser beams, this assumption is not false since the wings of a Gaussian beam profile extends, at least in theory, to infinity. Such an assumption would simplify the model M1. In the course of our studies, we assessed this approach, leading to a simplified model M3: with The results of a curve fitting with this simplified model M3 are shown in the graphs of Figure 16. The fitting was performed on the previous data using the limited range of the radial coordinate r as described above, see model M2. Although, the resulting fit parameters using model M3 are different to those of model M1, the course of the corresponding fitted curves is similar.
Sensors 2020, 20, x FOR PEER REVIEW 20 of 41 would simplify the model M1. In the course of our studies, we assessed this approach, leading to a simplified model M3: with ( , * ≔ 1) = 1 + .
The results of a curve fitting with this simplified model M3 are shown in the graphs of Figure  16. The fitting was performed on the previous data using the limited range of the radial coordinate as described above, see model M2. Although, the resulting fit parameters using model M3 are different to those of model M1, the course of the corresponding fitted curves is similar.

Summary of the Fit Procedure
For each measurement series, we performed the fit procedure as described before: The results of the fitting procedure are presented in detail in Appendix A and are summarized in the following Section 6.

Results
In this section, we present a summary of the results of the various lens scattering analyses and a subsequent statistical evaluation of the derived scatter parameters. Detailed results of the fitting procedure for each lens can be found tabulated in Appendix A. Figure 17 shows the results for the scatter parameters. In the left column, Figure 17a,c,e,g shows the values of the scatter parameters , , and as a function of the truncation factor , respectively. Each data point corresponds to a measurement series, as defined in Section 3. The data sets for the different lenses are distinguished by the color and shape of the data points. The numbering of the legend is equal to the numbering of the lenses as used in Tables 5 and 7. The scatter parameters , and were obtained directly from Fit 3, whereas scatter parameter was

Summary of the Fit Procedure
For each measurement series, we performed the fit procedure as described before: • Fit 1: Curve fitting with our original model M1 to the full pixel range of the measurement series. The results of the fitting procedure are presented in detail in Appendix A and are summarized in the following Section 6.

Results
In this section, we present a summary of the results of the various lens scattering analyses and a subsequent statistical evaluation of the derived scatter parameters. Detailed results of the fitting procedure for each lens can be found tabulated in Appendix A. Figure 17 shows the results for the scatter parameters. In the left column, Figure 17a,c,e,g shows the values of the scatter parameters s, b, l and b 0 as a function of the truncation factor ν, respectively. Each data point corresponds to a measurement series, as defined in Section 3. The data sets for the different lenses are distinguished by the color and shape of the data points. The numbering of the legend is equal to the numbering of the lenses as used in Tables 5 and 7. The scatter parameters s, b and l were obtained directly from Fit 3, whereas scatter parameter b 0 was calculated using the relation b 0 = b(100l) s ; see Equation (9). The scatter parameters are pure material parameters and are independent from the truncation factor ν. However, looking at these four graphs, we can see that the fluctuation of the values seems not always to be of statistical nature. In some cases, for example the scatter parameter l for lens #2 (Edmund Optics 67715), the data points seem to lie on a bended curve. We have no explanation for this observation, but attributed this to the simplicity of our theoretical models. The derivation of the model curve of Equation (1) was accompanied by several assumptions and simplifications to keep the equations manageable. A deeper investigation of this behavior would be of future interest. Furthermore, we can see some outliers in the data, for example for scatter parameter s regarding lens #5 (Navitar NMV-75). We therefore decided to use the median as a robust estimator for the central tendency of the scatter parameters for each lens. Table 7. Results of the statistical analysis for the scatter parameters s, b, l and b 0 . IQR: interquartile range, QCD: quartile coefficient of dispersion. The results of the statistical analysis for the scatter parameters s, b, l and b 0 are plotted as box plots in Figure 17b,d,f,h, respectively. The results are listed in Table 7: the median, the interquartile range (IQR, difference of third and first quartile) and additionally the quartile coefficient of dispersion (QCD). The QCD is defined as a ratio of IQR to the sum of the first and third quartile and describes the dispersion of the values. From the box plots of Figure 17 we can learn that for the camera lenses (lenses #1-#7), the medians of the scatter parameters s and b do not vary much. For the scatter parameter l, the variation of the median was somewhat stronger, caused mainly by lenses #2 (Edmund Optics 67715) and #7 (Schneider-Kreuznach Xenoplan 2.8/50). Additionally, for the scatter parameter b 0 , there was some larger fluctuation of the median value.  Figure 17. Results of the lens scattering analyses. The numbering of the data corresponds to the lens numbers given in Table 5 or Table 7. (a,c,e,g): Scatter parameters , , and as a function of the truncation factor . (b,d,f,h): Box plots of the scatter parameters , , and for the different lenses. Figure 17. Results of the lens scattering analyses. The numbering of the data corresponds to the lens numbers given in Table 5 or Table 7 The results of Table 7 build a promising basis to state a generic set of scatter parameters for COTS camera lenses with regard to analyze the incapacitation of sensors when dazzled. In other words, these results also help in terms of laser safety considerations or laser safety calculations as performed in reference [7]. Calculating the median of the scatter parameters of all camera lenses (lens #1-#7, excluding lens #8) results in the following generic set of scatter parameters for the data we measured: s = −1.88; b = 0.37 sr −1 ; b 0 = 6.74 sr −1 ; l = 2.02 mrad Unfortunately, these values do not fulfill Equation (9): b 0 = 6.74 sr −1 b(100 · l) s = 7.47. This relation is only valid for the scatter parameters obtained for each measurement series by the various curve fits (listed in Appendix A). Neither the median values for each data set (corresponding to a camera lens) stated in Table 7 do fulfill the relation nor does the generic set of scatter parameters stated above.
In order to be able to deduce a consistent set of scatter parameters, the values given above have to be adjusted accordingly. For example, one could simply increase scatter parameter s (to a less negative value), decrease scatter parameter b, increase b 0 or increase scatter parameter l until the relation of Equation (9) is valid, but this would be interpreted like a purely arbitrary approach. Thus, we chose another way by changing all scatter parameters slightly until the relation was fulfilled. Based on our own thoughts, the adaption of the scattering parameters was performed according to the equation where k is a constant. Equation (24) takes into account that the adaption of the scatter parameters by a constant k may have an influence larger than this factor. For example, a change of scatter parameter l by a factor of 0.98 would result in a change of a factor 0.96 on the right hand side of Equation (9) due to the exponent s. This is also true for scatter parameter s, which leads to prefactors of 1/ √ k and 3 √ k for parameters l and s, respectively.
Equation (24) is valid for a factor k ≈ 0.974, i.e., all parameters were changed by less than 2.6 percent. By rounding to the first two decimal places, we received a new generic set of scatter parameters, expressed now by the capital letters S, B, B 0 and L: S = −1.86; B = 0.36 sr −1 ; B 0 = 6.92; L = 2.04 mrad There may be alternative methods to balance the scatter values, but with Figures 18 and 19 we demonstrate the applicability of this method. The figures show the measured irradiance profiles for all seven camera lenses, plotted in separate graphs. Since the incident laser power P in varied for the different settings of the f-number F, we scaled all the profiles to the same (arbitrarily chosen) input power of 0.1 µW by applying the respective factor of 0.1 µW/P in to the data, in order to normalize the data. The resulting point cloud reflects quite well the value range of the camera signal µ y as a function of radial coordinate r px .
Furthermore, we plotted theoretical curves regarding the minimum and maximum f-number F min , F max of each lens. For the theoretical curves, we used the original model M1 (black curves) and additionally the simplified model M3 (blue dotted curve). The curve for model M3 was calculated only for f-number F min , since there was no difference to model M1 in the case of f-number F max . All theoretical curves were calculated using the reference wavelength of 550 nm. For the calculations, we used both the individual scatter parameters of the lenses as given by the median values of Table 7 (graphs on the left hand side) and the generic set of scatter parameters S, B and L as stated above (graphs on the right hand side).

Simulation of Stray Light Irradiance Using the Optical Engineering Software FRED
Our determination of a generic set of scatter parameters for camera lenses was based on measurements using a sample of seven different COTS camera lenses. It is clear that this set of values will not cover adequately all kinds of camera lenses. In order to compare the results of our theoretical stray light model (using the generic set of scatter parameters , , and ) with stray light distributions of other typical camera lenses, we utilized the optical engineering software FRED from Photon Engineering. Using this software, we performed stray light analyses for two camera lenses of We can see that the F min , F max model curves enclosed quite well the data points for radial coordinates r 10 pixel = 55 µm. It is clear that the curves calculated with the individual scatter parameters of the camera lenses gave more accurate results than the generic set of scatter parameters. However, the generic set seemed to be a good choice if scatter parameters of a camera lens are unknown. Depending on the camera lens, sometimes differences occurred between the result of model M1 and M3. In case of lenses no. 1, 2 and 7 (Edmund Optics 54690 and 67715, Schneider-Kreuznach Xenoplan 2.8/50), there was no difference visible. For lenses no. 3 and 6 (Edmund Optics 86410 and Navitar NMV-100), there appeared larger difference for radial coordinates r 50 pixel = 275 µm.
Regarding the sensor's incapacitation, this difference of models M1 and M3 may play a role for dazzle scenarios, especially when the laser beam diameter is smaller or similar to size of the lens aperture (ν < √ 2) and for rather small dazzle spots with r dazzle 50 pixels. For most practical applications, where the laser source is typically further away and the laser beam overspills the optics diameter, the simplified model M3 should be adequate.

Simulation of Stray Light Irradiance Using the Optical Engineering Software FRED
Our determination of a generic set of scatter parameters for camera lenses was based on measurements using a sample of seven different COTS camera lenses. It is clear that this set of values will not cover adequately all kinds of camera lenses. In order to compare the results of our theoretical stray light model (using the generic set of scatter parameters S, B, B 0 and L) with stray light distributions of other typical camera lenses, we utilized the optical engineering software FRED from Photon Engineering. Using this software, we performed stray light analyses for two camera lenses of the double Gauss type, since it is stated that "35-mm SLR normal lenses are invariably Double-Gauss types" [17]. Furthermore, we modeled the achromatic doublet lens (Thorlabs AC254-050-A) also used in our measurements.

Layout of the Stray Light Simulation
The lenses modeled with FRED are listed in Table 8 and the corresponding optical layouts are shown in Figure 20. The optical layouts for the double Gauss lenses were taken from reference [18]; the optical layout of the achromatic doublet lens was provided by the manufacturer. Table 8. Camera lenses modeled using the optical engineering software FRED. the double Gauss type, since it is stated that "35-mm SLR normal lenses are invariably Double-Gauss types" [17]. Furthermore, we modeled the achromatic doublet lens (Thorlabs AC254-050-A) also used in our measurements.

Layout of the Stray Light Simulation
The lenses modeled with FRED are listed in Table 8 and the corresponding optical layouts are shown in Figure 20. The optical layouts for the double Gauss lenses were taken from reference [18]; the optical layout of the achromatic doublet lens was provided by the manufacturer.  For the stray light simulation, we used the settings for the FRED software listed in Appendix C, Table A9. The laser beam was simulated by a source grid of a defined number of input rays, homogeneously distributed within a predefined aperture. The diameter of that aperture was adjusted in that way that the source grid was slightly larger than the maximum input aperture of the considered lens. The total input power of 2 µW was distributed to the rays of each source grid in such a way, to receive always a Gaussian beam with a -diameter of 21.1 mm. The wavelength of the simulated light source was set to 550 nm. These settings corresponded largely to those of our For the stray light simulation, we used the settings for the FRED software listed in Appendix C, Table A9. The laser beam was simulated by a source grid of a defined number of input rays, Sensors 2020, 20, 6308 27 of 42 homogeneously distributed within a predefined aperture. The diameter of that aperture was adjusted in that way that the source grid was slightly larger than the maximum input aperture of the considered lens. The total input power of 2 µW was distributed to the rays of each source grid in such a way, to receive always a Gaussian beam with a d 86 -diameter of 21.1 mm. The wavelength of the simulated light source was set to 550 nm. These settings corresponded largely to those of our experimental setup.
The simulation of light scattering comprised of two different kinds of scatter functions. First, scattering of light at the rough surfaces of the optical elements according to Harvey's scatter model and, second, scattering of light at the housing, at the rim of the apertures and at the lens' edges. For the second, we assumed Lambertian reflection of 4% reflectivity. As scatter parameters for the Harvey scatter model, we either used our generic set of scatter parameters as described above (lens Fr1/Fr2) or, in the case of the achromatic doublet (lens Fr3), the measured scatter parameters as stated in Table 7 (see results for lens #8). Since we had no exact CAD models of lens barrels for the double Gauss lens systems, their housing was simply simulated by conical tubes connecting the edges of adjacent optical elements. More sophisticated methods for stray light control, like baffles and vanes, special paints and surface treatments were not simulated. In the case of the achromatic doublet lens, a cylindrical tube was simulated as a lens barrel, which is close to reality since the housing for the achromatic doublet lens was built using a standard opto-mechanical tube system; see Figure 7h.
For the analysis of stray light in the focal plane of a sensor, we defined a detector of dimensions 5.5 mm × 5.5 mm with 1000 pixels × 1000 pixels. This results in a pixel size of 5.5 µm, which is identical to the pixel size of the imaging sensor we used in our measurements. The detector was always positioned at the geometric focus of the respective lens.
The basic principle of stray light simulations is that each of the input beams that will be refracted/reflected at the lenses' surfaces or mechanical parts according to geometric optics, generates a large number of scattered rays. Thus, the number of rays to be simulated in total is a high multiple of each single input ray. Consequently, the number of rays reaching the detector depends not only on the number of input rays, but also on the vast number of parameters of the simulation software, which influences the accuracy of the simulation. That is, not only the number of input rays but also the number of scattered rays have a major impact on the outcome. It can be said, that the larger the number of input rays, the better the results for the central spot that is dominated by diffraction and aberrations. On the other side, the scattered irradiance distribution can be simulated better when using a large number of scatter rays. Therefore, producing realistic results using such an optical engineering software requires both a huge number of input and scatter rays, which requires a huge amount of computational power. Here, we present first the results of our stray light simulation. Since the number of applicable rays was limited by the computer hardware, we decided to privilege the scattering part for this publication. This means that we expect and accept deviations between simulations and measurements in the regime where diffraction dominates the irradiance signal. More details on the simulation process and an extensive analysis of our results will be presented in a dedicated publication.

Simulation Results
In the course of our investigations, we observed that the simulations performed with or without lens housing nearly show the same results, as depicted in Figure 21 for the case of the achromatic doublet lens. The simulated irradiance data are plotted as a function of the radial distance with respect to the center of the laser spot for the two cases with lens housing (red data points) and without lens housing (blue data points). Both curves had more or less the same course, especially for larger values of the radial coordinate, where scattered light from the lenses itself dominated the signals in the focal plane. lens housing nearly show the same results, as depicted in Figure 21 for the case of the achromatic doublet lens. The simulated irradiance data are plotted as a function of the radial distance with respect to the center of the laser spot for the two cases with lens housing (red data points) and without lens housing (blue data points). Both curves had more or less the same course, especially for larger values of the radial coordinate, where scattered light from the lenses itself dominated the signals in the focal plane. This result can be attributed to the specific cases (paraxial setup) we investigated and we simulated here: light only impinges on the camera lens along the optical axis without hitting directly the housing. In turn this means that only scattered light from the optical elements will reach the This result can be attributed to the specific cases (paraxial setup) we investigated and we simulated here: light only impinges on the camera lens along the optical axis without hitting directly the housing. In turn this means that only scattered light from the optical elements will reach the housing and subsequently the sensor. In future, the effects of the oblique incidence of light shall be examined, where the housing is directly hit by the impinging light.
In Figure 22, we plotted the simulated irradiance data as a function of radial distance with respect to the center of the laser spot (colored data points) for the three modeled lenses. For comparison, the graphs show results for raytracing without any scattering (blue data points) and with lens scattering (green data points). All these data was simulated without housing since a full scattering investigation would be highly time consuming, and would not lead to a better understanding of details. Furthermore, we also plotted the result of our theoretical model M1 using the same parameters as we used for the FRED simulations. The black solid curve shows the irradiance E M1 (r) of model M1 according to Equation (20), which comprises diffraction and scattering of light. The black dashed curve solely shows the scatter part E s (r) according to Equation (11). The graphs on the left-hand side were results for the smallest f-numbers modeled; the graphs on the right-hand side show results for the largest f-numbers modeled.
For small values of the radial coordinate r, there was a larger difference between the simulated irradiance and the output of our theoretical model. This is unsurprisingly, since the simulation also comprised aberrations whereas our theoretical model did not. Furthermore, as stated above, we adjusted the simulation in that way that the results should be more exact for the scattered part of the simulated irradiance distribution.
For larger values of the radial coordinate, the simulated irradiance values were typically slightly below the values of the theoretical model. The result for the larger f-number of the 35 mm double Gauss lens in Figure 22d is an exception. In this case, the f-number was larger (F = 16) compared to the other lenses; the simulated irradiance values were slightly above the model curve. We do not exactly know where these deviations resulted from. This will be investigated in more detail in future work.
However, as a first result, the simulation is in good agreement with the theoretical model. Thus, we could conclude that the theoretical model applied with the generic set of scatter parameters allows for an appropriate estimate of the irradiance distribution of typical camera lenses.
with lens scattering (green data points). All these data was simulated without housing since a full scattering investigation would be highly time consuming, and would not lead to a better understanding of details. Furthermore, we also plotted the result of our theoretical model M1 using the same parameters as we used for the FRED simulations. The black solid curve shows the irradiance ( ) of model M1 according to Equation (20), which comprises diffraction and scattering of light.
The black dashed curve solely shows the scatter part ( ) according to Equation (11). The graphs on the left-hand side were results for the smallest f-numbers modeled; the graphs on the right-hand side show results for the largest f-numbers modeled.

Summary
In this publication, we present our measurements to assess the scattering parameters of commercial off-the-shelf (COTS) camera lenses. For this, the spatial irradiance distribution of laser light at the focal plane of different camera lenses was measured using a camera as a detector. Assuming rotational symmetry of the irradiance distribution, the image data was used to derive radial irradiance profiles. Subsequently, a simple theoretical model for irradiance distribution calculations was used to perform curve fitting to the experimentally measured irradiance profiles in order to extract the scatter parameters of the camera lenses. These scatter parameters are related to the well-known 3-parameter Harvey scatter model describing light scattering from the rough surfaces of optical elements. The main outcome of our work shows that the values of the scatter parameters for quite specific types of camera lenses were very similar. This allowed us to state a generic set of scatter parameters for typical COTS camera lenses and, moreover, now will allow us to perform laser safety calculations for sensors even in the case that exact values for the scatter parameters of a camera lens are not available, which was the main motivation for this work.
To our best knowledge, it is the first time that scatter parameters for standard COTS camera lenses were published. However, the model used to extract the scatter parameters was not a rigorous theoretical model of the irradiance distribution within camera lenses. The model was specifically developed to perform laser safety calculations for imaging sensors. It comprised several simplifications and, thus, the scatter parameters presented here are specific and might not be readily applicable for other types of optics. Therefore, we used the FRED optical engineering software to examine whether the stated general set of scatter parameters can also be applied to other camera lenses than those used for the measurements presented here. As a result, we see that in combination with the dedicated theoretical model, the stated general set of scatter parameters allows a good estimation of the focal plane irradiance distribution of camera lenses.
Future work on that topic could comprise the increase of the statistical database by both including further camera lenses and performing more measurements on the currently used camera lenses. Regarding further camera lenses, additional investigations could include, for example, zoom lenses, where the scattering parameters could be measured at different settings of their focal length, or a telephoto lens with a very large focal length (500 mm). For the camera lenses already in use, further research could comprise measurements for different truncation factors by changing the laser beam diameter instead of the f-number. It would also be of interest to test how non-coherent radiation (e.g., by using light of a narrowband LED) would change the stray light distribution. An important point for future work is the validation of laser safety calculations for imaging sensors using the here stated generic set of scatter parameters. The validation could be performed using data acquired with different sensors and different camera lenses, in particular data gathered in free-field conditions. Furthermore, it is of particular interest whether the simplified model already meets the requirements for the laser safety calculations. This would be another step to a further simplification of our model.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Results of the Data Analysis
As a result of our data exploitation, Tables A1-A8 present the scatter parameters s, b and l, which served as curve fitting parameters to fit our data by means of equations derived from our models M1-M3. Additionally, the scatter parameter b 0 , calculated from Equation (9), is presented. The fit procedure, explained in detail in Section 5.3, was based on four individual curve fittings. For the convenience of the reader, the various curve fits are repeated here: The measurements were performed on eight different lenses and for various settings of f-number F/mean truncation factor ν. All results were converted to the reference wavelength of 550 nm. Each table also contains the median of the various scatter parameters. The values of fit no. 3 are shown in bold letters, since these values were used for the statistical analysis of Section 6. To calculate the median values for a specific lens, usually all the fit results for this lens were taken into account. However, in some cases the fit value for the scatter parameter l reached its lower limit of 0.5 mrad. These cases are highlighted by the red letters in the tables. When this occurred for fit no. 3, the results for this measurement series were excluded from the calculation of the median values, since these fit results were assumed to describe the data not properly. The excluded data is highlighted by the orange background in the tables. The determination of the laser beam diameter at the entrance aperture of the camera lenses was an important task that also required great care. Since commercially available beam profilers are limited to smaller laser beam sizes, we measured the beam diameter using the setup depicted in Figure A1. The setup was largely equal to the setup of Figure 4, except that the camera C and camera lens CL were replaced by a white, diffuse viewing screen VS (Thorlabs EDU-VS1/M). Furthermore, a beam-monitoring camera BMC (camera Allied Vision Mako G-223B NIR + camera lens Edmund Optics 67716) was placed near the off-axis parabolic mirror OPM to observe the viewing screen.
Before the measurement of the laser beam diameter, we ensured that the BMC was focused to the viewing screen. Furthermore, we estimated the magnification M bmc for the setup by placing a resolution test chart at the position of the viewing screen. Using the known size of the test chart's features and the corresponding size in the camera image, the magnification was calibrated to be M bmc = 124.6 µm/pixel. The determination of the laser beam diameter at the entrance aperture of the camera lenses was an important task that also required great care. Since commercially available beam profilers are limited to smaller laser beam sizes, we measured the beam diameter using the setup depicted in Figure A1. The setup was largely equal to the setup of Figure 4, except that the camera C and camera lens CL were replaced by a white, diffuse viewing screen VS (Thorlabs EDU-VS1/M). Furthermore, a beam-monitoring camera BMC (camera Allied Vision Mako G-223B NIR + camera lens Edmund Optics 67716) was placed near the off-axis parabolic mirror OPM to observe the viewing screen.
Before the measurement of the laser beam diameter, we ensured that the BMC was focused to the viewing screen. Furthermore, we estimated the magnification for the setup by placing a resolution test chart at the position of the viewing screen. Using the known size of the test chart's features and the corresponding size in the camera image, the magnification was calibrated to be = 124.6 µm/pixel. We captured camera images of the laser beam for all four laser wavelengths and dark frames, i.e., without laser illumination. In each case, 10 camera images were acquired and averaged before further analysis. Figure A2 shows the mean images of the laser beam profiles for all four wavelengths. For each image, the centroid of the irradiance distribution was estimated and then a vertical and horizontal profile at the centroid's position was extracted. For each vertical/horizontal profile, 21 pixel columns/rows were used and averaged. These areas are marked by the red bars in the images of Figure A2. The corresponding profile lines are plotted in white color into the images. By fitting Gaussian curves to the profile lines, the laser beam diameter for each wavelength was estimated. The values for the vertical and horizontal profile and the mean values are printed in each image. We captured camera images of the laser beam for all four laser wavelengths and dark frames, i.e., without laser illumination. In each case, 10 camera images were acquired and averaged before further analysis. Figure A2 shows the mean images of the laser beam profiles for all four wavelengths. For each image, the centroid of the irradiance distribution was estimated and then a vertical and horizontal profile at the centroid's position was extracted. For each vertical/horizontal profile, 21 pixel columns/rows were used and averaged. These areas are marked by the red bars in the images of Figure A2. The corresponding profile lines are plotted in white color into the images. By fitting Gaussian curves to the profile lines, the laser beam diameter d 86 for each wavelength was estimated. The values for the vertical and horizontal profile and the mean values are printed in each image. Using the magnification M bmc , the mean laser beam diameter d 86 was estimated to be 21.53 mm, 21.12 mm, 20.30 mm and 20.14 mm for the laser wavelength 488 nm, 515 nm, 561 nm and 640 nm, respectively. Since the beam profiles did not exhibit a distinctive ellipticity, we used these mean beam diameters for the data analysis of our stray light measurements (see Section 5). Using the magnification , the mean laser beam diameter ̅ was estimated to be 21.53 mm, 21.12 mm, 20.30 mm and 20.14 mm for the laser wavelength 488 nm, 515 nm, 561 nm and 640 nm, respectively. Since the beam profiles did not exhibit a distinctive ellipticity, we used these mean beam diameters for the data analysis of our stray light measurements (see Section 5).

Appendix B.2. Characterization of the Camera
For the analyses of the images acquired during the experiments, it was essential to characterize the camera before the measurements. This comprised the verification of linearity and the estimation of its characteristics (see Table 2): saturation gray value . , absolute sensitivity threshold . , overall system gain and the quantum efficiency . The characterization was conducted in accordance with the EMVA 1288 standard [15]. However, we performed only measurements necessary to estimate the aforementioned parameters and not the complete procedure for a valid camera characterization corresponding to the EMVA 1288 standard. Figure A3 shows a schematic drawing of the experimental setup for the characterization. The illumination of the imaging sensor was provided by a light-emitting diode LED (Thorlabs M625F2) connected to an integrating sphere IS (Labsphere 3P-GPS-033-SL). In accordance with EMVA 1288 standard, the camera C was located at a distance to the exit port of the integrating sphere that equals eight times the exit port's diameter. To prevent illumination of the imaging sensor from extraneous light, the beam path between the exit port of the IS and the camera was enclosed by a tube. Additionally, the lighting of the laboratory was switched off during measurements. A second exit port of the IS was used to connect a bifurcated fiber bundle (Thorlabs BFY400LS02). The two legs of the fiber bundle were connected to a reference photodiode PD (Ophir PD300R-UV) and a spectrometer SM (RGB Photonics Qwave).

Appendix B.2. Characterization of the Camera
For the analyses of the images acquired during the experiments, it was essential to characterize the camera before the measurements. This comprised the verification of linearity and the estimation of its characteristics (see Table 2): saturation gray value µ y.sat , absolute sensitivity threshold µ y.min , overall system gain K and the quantum efficiency η. The characterization was conducted in accordance with the EMVA 1288 standard [15]. However, we performed only measurements necessary to estimate the aforementioned parameters and not the complete procedure for a valid camera characterization corresponding to the EMVA 1288 standard. Figure A3 shows a schematic drawing of the experimental setup for the characterization. The illumination of the imaging sensor was provided by a light-emitting diode LED (Thorlabs M625F2) connected to an integrating sphere IS (Labsphere 3P-GPS-033-SL). In accordance with EMVA 1288 standard, the camera C was located at a distance to the exit port of the integrating sphere that equals eight times the exit port's diameter. To prevent illumination of the imaging sensor from extraneous light, the beam path between the exit port of the IS and the camera was enclosed by a tube. Additionally, the lighting of the laboratory was switched off during measurements. A second exit port of the IS was used to connect a bifurcated fiber bundle (Thorlabs BFY400LS02). The two legs of the fiber bundle were connected to a reference photodiode PD (Ophir PD300R-UV) and a spectrometer SM (RGB Photonics Qwave). Before the actual measurements were performed, the reference photodiode PD had to be calibrated by placing a power meter at the position of the camera's imaging sensor and measuring the ratio of the signals when the LED was switched on. Using this ratio, we were able to calculate the irradiance at the position of the imaging sensor during the measurement by the readings of the reference photodiode. Furthermore, we could monitor whether the irradiance was kept constant during the measurement. Using the spectrometer SM, we measured a peak emission wavelength of 635 nm and a full width at half maximum (FWHM) spectral width of 20 nm for the LED.
For the characterization of the camera, we acquired a series of images for 56 different settings of the camera's exposure time from 2800 µs (all pixel saturated) to the minimum of 60 µs in steps of 50 µs. For each setting, four images were acquired: two each with and without irradiation (dark frames). Using this image data, we calculated the mean gray value and the temporal variance of the gray values for the illuminated images and the corresponding quantities . and . for the dark frames. The specification to calculate these quantities is given by Equations (28) and (29) of the EMVA 1288 standard [15]. Furthermore, the mean number of photons per pixel was calculated using Equation (2) of the EMVA 1288 standard: = /(ℎ / ). Using these quantities, the necessary camera parameters could be derived. Figure A4 shows a plot of the so-called photon transfer curve, where the photo-induced variance − .
is plotted as blue data points versus the mean photo-induced gray value − . . From this plot, the saturation gray value . and the overall system gain of the camera can be extracted. The saturation gray value is given as the mean gray value where the variance − .
has its maximum (marked by a green data point in the plot); which is . = 3861 DN for our camera. The overall system gain is given by the slope of the curve. We estimated the overall system gain to be = 0.399261 DN/e by a fit of a linear curve to the photon transfer curve (dashed line in the plot). The data range used for the fit was the data between the minimum value of − . and 0.7 ⋅ . , marked by red data points in the plot. In Figure A5, the sensitivity curve of the camera is shown. In this figure, the mean photo-induced gray value − . is plotted versus the mean number of photons arriving at a pixel. The slope of the linear part of this curve is the responsivity of the camera. This value can also estimated by the fit of a linear curve to the data using the same data range as before and was estimated to be = 0.30595 DN/photon. The responsivity and the overall system gain are related by = . Using this relation the quantum efficiency at the wavelength of 635 nm was calculated to be ≈ 0.77 e /photon.

Using the number of photons
corresponding to the saturation gray value . , we could calculate the saturation capacity  Before the actual measurements were performed, the reference photodiode PD had to be calibrated by placing a power meter at the position of the camera's imaging sensor and measuring the ratio of the signals when the LED was switched on. Using this ratio, we were able to calculate the irradiance at the position of the imaging sensor during the measurement by the readings of the reference photodiode. Furthermore, we could monitor whether the irradiance was kept constant during the measurement. Using the spectrometer SM, we measured a peak emission wavelength of 635 nm and a full width at half maximum (FWHM) spectral width of 20 nm for the LED.
For the characterization of the camera, we acquired a series of images for 56 different settings of the camera's exposure time from 2800 µs (all pixel saturated) to the minimum of 60 µs in steps of 50 µs. For each setting, four images were acquired: two each with and without irradiation (dark frames). Using this image data, we calculated the mean gray value µ y and the temporal variance of the gray values σ 2 y for the illuminated images and the corresponding quantities µ y.dark and σ 2 y.dark for the dark frames. The specification to calculate these quantities is given by Equations (28) and (29) of the EMVA 1288 standard [15]. Furthermore, the mean number of photons µ p per pixel was calculated using Equation (2) of the EMVA 1288 standard: µ p = AEt exp /(hc/λ). Using these quantities, the necessary camera parameters could be derived. Figure A4 shows a plot of the so-called photon transfer curve, where the photo-induced variance σ 2 y − σ 2 y.dark is plotted as blue data points versus the mean photo-induced gray value µ y − µ y.dark . From this plot, the saturation gray value µ y.sat and the overall system gain K of the camera can be extracted. The saturation gray value is given as the mean gray value where the variance σ 2 y − σ 2 y.dark has its maximum (marked by a green data point in the plot); which is µ y.sat = 3861 DN for our camera.
The overall system gain is given by the slope of the curve. We estimated the overall system gain to be K = 0.399261 DN/e − by a fit of a linear curve to the photon transfer curve (dashed line in the plot). The data range used for the fit was the data between the minimum value of µ y − µ y.dark and 0.7 · µ y.sat , marked by red data points in the plot. In Figure A5, the sensitivity curve of the camera is shown. In this figure, the mean photo-induced gray value µ y − µ y.dark is plotted versus the mean number of photons µ p arriving at a pixel. The slope of the linear part of this curve is the responsivity R of the camera. This value can also estimated by the fit of a linear curve to the data using the same data range as before and was estimated to be R = 0.30595 DN/photon. The responsivity R and the overall system gain K are related by R = ηK. Using this relation the quantum efficiency at the wavelength of 635 nm was calculated to be η ≈ 0.77 e − /photon. Using the number of photons µ p corresponding to the saturation gray value µ y.sat , we could calculate the saturation capacity µ p.sat ≈ 12, 488 photons. The saturation capacity of our camera related to electrons is µ e.sat = η · µ p.sat ≈ 9570 e − .
From the sensitivity plot of Figure A5, we could also see that the linearity of the camera is quite good, which means that the camera was well suited for the kind of measurements presented in this publication.
Sensors 2020, 20, x FOR PEER REVIEW 39 of 41 From the sensitivity plot of Figure A5, we could also see that the linearity of the camera is quite good, which means that the camera was well suited for the kind of measurements presented in this publication.    Table 7, results for lens #8) From the sensitivity plot of Figure A5, we could also see that the linearity of the camera is quite good, which means that the camera was well suited for the kind of measurements presented in this publication.    Table 7, results for lens #8)