The Development of a 3D LADAR Simulator Based on a Fast Target Impulse Response Generation Approach

A new laser detection and ranging (LADAR) simulator has been developed, using MATLAB and its graphical user interface, to simulate direct detection time of ﬂight LADAR systems, and to produce 3D simulated scanning images under a wide variety of conditions. This simulator models each stage from the laser source to data generation and can be considered as an efﬁcient simulation tool to use when developing LADAR systems and their data processing algorithms. The novel approach proposed for this simulator is to generate the actual target impulse response. This approach is fast and able to deal with high scanning requirements without losing the ﬁdelity that accompanies increments in speed. This leads to a more efﬁcient LADAR simulator and opens up the possibility for simulating LADAR beam propagation more accurately by using a large number of laser footprint samples. The approach is to select only the parts of the target that lie in the laser beam angular ﬁeld by mathematically deriving the required equations and calculating the target angular ranges. The performance of the new simulator has been evaluated under different scanning conditions, the results showing signiﬁcant increments in processing speeds in comparison to conventional approaches, which are also used in this study as a point of comparison for the results. The results also show the simulator’s ability to simulate phenomena related to the scanning process, for example, type of noise, scanning resolution and laser beam width.


Introduction
Laser detection and ranging, or laser radar (LADAR) systems, are considered an attractive alternative to radio detection and ranging (RADAR) systems because they use laser wavelengths which are shorter than RADAR wavelengths, to produces very highresolution 3D images. In addition, light velocity allows LADAR systems to take numerous measurements per second. LADAR images are created by scanning a scene with laser beams, the return time for these beams used to calculate range LADAR data. The format for this data is range, azimuth and elevation angle, this representing the spherical coordinates system whose origin is the sensor. LADAR converts this type of data into a 3D Cartesian format in order to produce a three-dimensional range image, this in turn representing the spatial location of the intersection of the laser beam with the scanned scene.
LADAR systems play diverse roles in both civilian and military applications. In ground navigation, they are used for obstacle and road-boundary detection, and autonomous vehicle navigation [12,14,19,29,31,32,44,47]. In aerial navigation, they provide autonomous navigational capacities [16], obstacle warning systems [11], and considered as a reliable alternative to GPS [40]. Regarding maritime navigation, LADAR are used for both precise manoeuvring operations and obstacle avoidance [24]. Looking to their use by the military, LADAR assists target detection and classification [8,10,39], anti-ship missile tracking [33], target identification at long range [5-7, 25, 41, 42], and the identification of military ground vehicles that may be hidden under camouflage or foliage such as tree canopies [30].
In consequence, simulations for LADAR systems have become a valuable tool for developing said systems and their data processing algorithms [13] because the simulator is able to produced LADAR images under different controlled effects, which enable the algorithms' developers to evaluate their algorithms under these effects individually. In order to simulate these systems, a target impulse response must be generated for each laser pulse transmitted to the components of the target. This is an extensive computational process as the intersection points of each laser beam with the target's surface, need to be identified in addition to their traveling distances.
Since most applications related to developing LADAR systems require LADAR simulators to be able to accurately simulate the propagation of the laser beam very quickly; this implies the need for rapid processing of a large number of laser footprint samples. In addition, in order to develop LADAR processing algorithms using simulated data, the simulators must be able to scan a large number of targets at high speed, under different scanning parameters.
Several methods have been developed to increase the speed of the required computational process; some approximate the simulation by defining the reflection of the laser pulse as 3D model voxels that have a direct lineof-sight to the sensor [20,21]. Others calculate the distance between the target and the sensor by using the division of the model's surface and the distance between the viewpoint and model's 3D points for simplification [45,46]. A single wide laser beam projection, with focalplane array and parallel computing, is also used to reduce computational time [15,22,26,43].
In this paper, a new approach to generate the actual target impulse response, based on finding the actual intersection points between the laser beam and the 3D model, is presented. This approach is based on deriving the algorithms required to calculate target angular ranges, these algorithms used to speed up the process. In order to evaluate the performance of the simulator using this approach, over forty 3D models were scanned, under different scanning parameters, the simulation times recorded.
In the following sections, the theoretical background, which includes the equations and parameters required to simulate the laser beam propagation and thus the core of the simulator, are presented. This is followed by the new approach of generating the target impulse response. The main concepts of the simulation implementation for the LADAR simulator and its control windows are described using a collection of selected simulated images. The testing procedure and results of the evaluation are given followed by the conclusion.

Simulation of Laser Beam Propagation
The process of simulating the propagation of the laser beam from the transmitter to the target and back to the receiver is presented in this section in order to provide a clear understanding of the target impulse response (TIR). This process is divided into four parts, as shown in Fig. 1 [38]. The laser beam energy distribution is described first in both temporal and spatial domains

Laser Beam Energy Distribution
In order to simulate laser beam energy distribution, the laser pulse has to be modelled in time (temporal distribution) as well as in space (spatial distribution) resulting in a four-dimensional model. The outgoing pulse intensity is decomposed as [34]: Gðt; H ls ; V ls ; R ls Þ ¼ pðtÞ Â IðH ls ; V ls ; R ls Þ ð 1Þ where p(t) is the discrete pulse shape in time domain; IðH ls ; V ls ; R ls Þ is the proportion of energy contained within a component at a location of H ls ; V ls ; R ls dimensions, where t; H ls ; V ls ; R ls take on discrete values. Figure 1 shows these two distributions for the laser pulse travelling from the LADAR system to the target, where H ls and V ls are the horizontal and vertical cross-range dimensions respectively and R ls is range dimension in the direction of the pulse travelling.

Temporal Distribution
The amount of laser power that is produced by the LADAR source and transmitted toward the target area is assumed to have a Gaussian distribution with time. This distribution is defined by the laser pulse energy E t in a unit of joules and pulse width s (full-width at halfmax power in a unit of seconds). The following equation [13,34] is used to model this distribution: where p(t) is the laser pulse power in unit of watts at time t. In order to avoid aliasing on the laser pulse waveform, the sampling period Dt is calculated using the following equation [34]: where sf t is the time sampling factor. Its values must be greater than one for good pulse representation.

Spatial Distribution
The laser intensity profile produced by a laser source cavity is not constant across the beam diameter at all ranges, and it depends on the technique used to generate the laser beam. Generally, this profile or energy distribution is modelled as a spatial Gaussian function in horizontal H ls , vertical V ls , and range R ls dimensions [27,36]: where WðR ls Þ is the beam width (radius) at R ls defined as the radial distance (in metres) at which the profile value is decreased to 1=e 2 from its peak value. The dependence of beam width on the distance R ls is governed by [27,36]: where W o is the beam waist at R ls ¼ 0 and k is the laser wavelength.
Referring to [36] and [27], the spatial sampling size Ds that avoids aliasing of the laser intensity profile, can be calculated using the following equation: where sf s is the spatial sampling factor. Its values must be greater than one for good intensity representation.

Atmospheric Effects
When the transmitted laser beam propagates through the atmosphere, some of the energy is absorbed and scattered by atmospheric molecules, dust and aerosols [34].This attenuation limits the performance of the LADAR system and is dependent on the laser's wavelength k and propagation length R ls . This can be modelled using Beer's law as shown by the following equation [23,34]: where T a is the one-way atmospheric transmission value and r s ðkÞ is the atmospheric coefficient in m À1 for the wavelength k.

Target Interaction
The interaction between the transmitted laser beam and the target surface produces a reflected signal. The characteristics of this signal depend on the surface reflectance q tr (2 À 25%) [34], the angle of dispersion X tr (Lambertian targets are assumed i.e. X tr ¼ p [9,23,34]), surface area A tr (extended targets are assumed), surface shape and beam incidence angle. The beam-target interaction not only affects the energy in the reflected signal, but also impacts on its shape. Both the beam incidence angle and surface shape (plane and step etc.) effects, play important roles in changing the transmitted pulse shape. To simulate these effects, the target impulse response (TIR) must be generated. The simulated shapes for the reflected signals that result from the interaction of the laser beam at different incident angles with the step targets are shown in Fig. 2.

LADAR Receiver
The process of determining the range to the target from the reflected signal is accomplished by the LADAR receiver. This process depends on detection techniques [3] (direct or coherent), optical transmission T o (the fraction of energy that arrives at the detector from the total energy captured by the receiver aperture), quantum efficiency g (the fraction of the signal that is converted into photoelectrons) of the detector, pulse detection technique and receiver noise (photon counting, speckle noise and background noise) (see ''Appendix 1'').
In this simulator, the direct detection technique is used by which the received optical energy is focused onto a photodetector element, with constant fraction discrimination (CFD) pulse peak detector, as this is unaffected by the amplitude fluctuation that causes jitter at the time of arrival [3]. The received signal power at the receiver aperture P r can be calculated using modified LADAR range equation (see ''Appendix 2'') [3,17,34]: where P t is the transmitter pulse power and D r is the diameter of circular receiver aperture.
More advanced and complex models [13,17,35,38,48,49], can be used to simulate the propagation of the laser beam, for example, from laser energy distribution to atmospheric effects and beam interaction models, to receiver optics and received signal processing electronics models. However, this study is focusses on evaluating the processing speed of the Fig. 2 Illustration of return pulse shaping by step target geometries [4] proposed TIR approach when used with the LADAR simulator. Standard laser propagation models are used as this will preserve the generality and give a clear indication about performance under fundamental (standard) models.
3 Proposed Approach to Generating The Target Impulse Response In order to generate a target impulse response for every laser burst, a laser beam footprint that illuminates the target surface IðH ls ; V ls ; R ls Þ must be created by using Eq. 4. The reflected power (P sample i ) reaching the receiver from each sample in this laser footprint and the corresponding round-trip time (i), are then calculated.
The sample reflected power is calculated from the LADAR range Eq. 8, while the round-trip time is calculated by obtaining the intersection for this sample with the target's surface. The reflected sample powers P sample i are then summed with the same time indices i to create the target impulse response (h tr ).
To obtain the intersection points of ray vectors representing the laser beam samples with triangular faces (see Fig. 3) that represent the model surface, each point (P) can be defined by the following equation [28]: where S represents the ray's starting position, V represents the direction in which the ray points, N is the triangle plane normal and P 0 , P 1 , and, P 2 are the triangle's vertices. If the denominator is equal to zero, then no intersection occurs. The barycentric coordinates [18] for this point, with respect to the triangle's vertices, are then calculated to determine if it lies inside the triangle's edges. In general, a point is inside (or on) the triangle if, and only if 0 w 1 1; 0 w 2 1, and w 1 þ w 2 1.
By defining the following vectors v 0 ¼ P 1 À P 0 ; v 1 ¼ P 2 À P 0 ; v 2 ¼ P À P 0 and using Cramer's rule, the w 1 , and w 2 values can be obtain through the following equations: To get the total laser footprint samples on the target model, the procedure of calculating the intersection point must be applied between every laser ray vector and all of the model's triangles [37]. This conventional (normal) approach is simple and straightforward but the number of calculations required make it a very time consuming approach especially when the model consists of a large number of triangles, scanning with high resolution, or when a large number of laser footprint samples are required.
To overcome these limitations, another approach is proposed. This approach is to only evaluate the intersection points between the vectors and triangles that lie in the same angular range; this allows a reduction in the number of calculations required thus speeding up the process. This novel approach preserves the fidelity that accompanies the increment in speed and produces a TIR identical to that generated from the previous approach, but by using less computational time. Figure 3 shows the principles used, where the steps in the procedure are presented as follows: 1. The angular extant in terms of azimuth and elevation angular ranges for each triangle is calculated and stored. These calculation are required ones per scanning setup. 2. Laser ray vectors (right side of Fig. 3) are generated. These vectors depend on the LADAR viewing direction, laser footprint size and the number of laser footprint samples. 3. The triangles whose angular extents (calculated in step 1) lie within the laser beam illumination direction, are selected (the blue edges triangles in Fig. 3). 4. For every selected triangle, the laser ray vectors that lie in the field of that triangle are selected and the intersection points between each calculated using Eqs. 9, 10, and 11. The top right side of Fig. 3, shows the selected rays that lie in the field of the green edged triangle. It also shows the intersection points on the triangle plane (green & yellow points) and inside the triangle itself (green points). 5. If the laser ray vector lies in the field of more than one triangle and has intersection points with each, the point that has the shorter distance to the laser is selected and stored.
The azimuth and elevation angular ranges mentioned previously (step 1) are calculated as follows: • Azimuth angular range: This is computed by calculating the azimuth angle for each triangle's vertices and comparing these angles with each other to find the minimum and the maximum values, these representing the azimuth angular range. • Elevation angular range: The method for calculating this range is similar to the method above except that the elevation angles for the triangle's vertices do not always represent the range. Therefore, additional three edge angles (one per triangle edge) are calculated and added to the comparison. Figure 4a shows the calculation principle for the edge angle / l . It starts by finding the line equation for the triangle edge of points P 0 & P 1 (red line in Fig. 4a) by: where P l is any point in the line of parameter l l . Its elevation angle / l can be calculated by: The first derivative for Eq. 13 is then taken and solved for zero, in order to find the parameter value l rp at which there is a round point P rp . After the derivation and simplification of the Eq. 13 it becomes: where U1 ¼ ðz 0 x 2 1 Þ þ ðz 0 y 2 1 Þ þ ðz 1 x 2 0 Þ þ ðz 1 y 2 0 Þ U2 ¼ Àðz 0 x 0 x 1 Þ À ðz 0 y 0 y 1 Þ À ðz 1 x 0 x 1 Þ À ðz 1 y 0 y 1 Þ If l rp is between 0 and 1, an additional elevation angle is required. Its value can be calculated by substituting l rp value into Eq. 13. Figure 4b shows the vertices elevation angles (at the start and at the end of curve) and the additional edge angle at the round point P rp (middle red circle).

LADAR Simulator
In this software, the simulated LADAR image is produced by scanning the model with multiple laser pulses, each providing one pixel on the LADAR image. The simulation steps for this simulator are as follows: 1. The required simulation parameters are defined; LADAR (viewing direction, field of view and scanning resolution), laser source (temporal and spatial domains), atmosphere, target, noise and receiver. 2. For every laser pulse, the target impulse response h tr is generated and convolved with the temporal laser pulse p(t) to calculate the temporal reflected power signal arriving at the detector P r ðtÞ as shown in the following equation.
3. The resultant power signals are converted to photoelectrons using Eq. 18. The background, photon counting, and speckle noise are then applied, if enabled by the user, using Eqs. 19 and 17 (see ''Appendix 1''). 4. The received electrical signals are then passed to the CFD peak detector to detect their peaks which are then used to calculate the round-trip time intervals Dt tot as shown in Fig. 1. As the laser pulses travel at the speed of light c (3 Â 10 8 m / s), the LADAR system calculates the ranges R c for these pulses using the following equation [34]: 5. Finally the range values are assigned to the corresponding pixels on the LADAR image.
Two graphical user interface (GUI) Windows were designed for this simulator. The main window (Fig. 5), is used to set the scan parameters and display the results. The second window (Fig. 6), is used to adjust the position and orientation for both the LADAR and the model (model size also can be adjusted in this window). The simulator is designed to display the results after the scanning process in four different formats, one in spherical coordinates (spherical image), the other

Testing Procedure and Evaluation Results
The performance of the LADAR simulator, in terms of processing time required to generate simulated 3D images, has been evaluated. This evaluation was achieved with the simulator using two target impulse response TIR generation approaches; conventional (normal) and proposed. Since the time required to generate TIR depends on the number of rays' vectors (that represent the laser beam samples) and on the number of triangles (that represent the scene or model surface), the simulator only tested changes in the effects of these parameters. The other scanning parameters were kept constant, their values presented in Table 1 in ''Appendix 3''. This table also shows the specifications for the computer that is used to run this test. Changing the number of vectors is achieved by changing the spatial sampling factor (sf s , see Sect. 2.1.2), while changing the number of triangles is done by using different 3D models of different numbers of faces.
To cover a wider range, more than forty 3D models of triangles with numbers ranging from 2375 to 13980, were used. Some models were taken from 3DVIA and ARTIST 3D model libraries [1,2] (see Fig. 9), the others generated from scanning real objects with triangulation based LADAR systems (TriLADAR). Figure 10 shows two TriLADAR systems with their control windows. These systems were designed and implemented at Liverpool University, to scan real objects and generate their 3D models with a specific number of triangles.
The testing procedure starts by scanning the 3D model using both approaches at a scanning resolution equal to 2500 pixels, with a laser beam (number of vectors equal to 768 by setting sf s to 5), and calculates the required time to get the final image. The procedure then increases sf s by 5 and re-scans the model again until the sf s reaches 50 this equivalent to 68403 vectors. Another 3D model comprised of more triangles than the previous model is then selected and the whole procedure is repeated again and so on, until 42 different 3D models are scanned. In order to guarantee that all models are fully scanned with the same resolution (2500 pixels), both the angular field of view and angular resolution are automatically adjusted according to the model dimensions.
The full testing results representing the effects of changing the number of triangles and vectors on the execution time for both approaches (normal and proposed), are presented as a 3D-graph as shown in Fig. 11a. Since the execution time for the normal approach covers a wider range compared to the proposed approach, the execution time is plotted in logarithmic scale as shown in Fig. 11b.
In order to present these effects individually, some results have been selected from the original 3D-graph and their slopes are also calculated (using the least square method) as shown in Fig. 11c, d. Figure 11c shows the effect of changing the number of triangles on execution time (and its slope) for specific vectors numbers (Vc. No.) while Fig. 11d shows the effect of changing the number of vectors on execution time (and its slope) for specific numbers of triangles (Tr. No.).   The results in Fig. 11b shows that the execution time for the proposed approach is much smaller than the normal approach. Figure 11c shows an increment in slopes for both approaches when the numbers of vectors increases from 768 to 68403. The differences in execution time when increasing the numbers of triangles, are not significant with the proposed approach. This caused the effect of models shapes clearly noticeable as a fluctuation in time (see Fig. 11c). Figure 11d also shows an increment in slopes for both approaches, but this time when the numbers of triangles increases from 2375 to 13980. In general, the average execution times for both normal and proposed approaches, are equal to 6:1 Â 10 3 s and 17.6 s respectively.

Conclusion
An efficient LADAR simulator has been developed using a novel TIR generation approach, to simulate the direct detection time of flight LADAR systems. The simulator models each stage, from laser source to data generation, over a short execution time producing simulated LADAR images, under a wide variety of conditions. The proposed approach to generate TIR has been developed to produce responses identical to these generated from the conventional or standard approach, but by using less computational time. This has been achieved by mathematically deriving the required equations to calculate target angular ranges which, in turn, enables an evaluation of the intersection points that lie in the same angular range, instead of evaluating the whole intersection point (between every laser ray vector and all the scene's triangles).
More than forty, 3D models were used to evaluate the simulator's performance in terms of processing time with different laser beam samples. The evaluation was carried out with this simulator using two target impulse response TIR generation approaches, proposed and conventional (normal), where the latter was used to benchmark the results. A comparison of the results shows that the LADAR simulator with proposed approach is quicker than normal approach, especially when the 3D model consists of a large number of triangles, or when a large number of laser footprint samples are required.
The average processing speed for the simulator with the proposed approach was 345 times faster in comparison to the normal one. This improvement in performance enables the simulator to scan a large number of targets, at different scanning parameters and poses at high speed, and opens up the possibility for simulating LADAR beam propagation more accurately in a shorter time, by using a large number of laser footprint samples.
The simulation steps for the LADAR simulator and its GUI are illustrated with some results of scanned 3D models. These simulation results demonstrate the ability of the LADAR simulator to scan and produces LADAR images under different scanning parameters (noise type, scanning resolution and laser beam width).
Acknowledgements I would like to express my sincere gratitude to those who have provided the financial support throughout the research, the Iraqi Ministry of Higher Education and Scientific Research. Also I would like to thank Dr. Toby Hall (Mathematical Sciences Department in University of Liverpool) for his help.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix 1: Photon Counting, Speckle, and Background Noise
The photon counting noise is defined as the random arrival of the reflected photons to the detector. This randomness introduces uncertainty in the number of photons measured during a finite time interval Dt. While the laser speckle noise is caused by the electromagnetic field interference occurring at the detector surface from a large collection of independent coherent radiators.
The statistical fluctuations in the measured signal due to speckle and photo counting noise can, be simulated in the LADAR measurement by modelling the number of photons detected as a negative binomial random variable with a mean equal to E½N signal and a variance given by the following equation [34]: where r 2 noise is the variance of the measured photocounts, and M CH is the degrees of freedom number for the laser light. For fully coherent light, M CH ¼ 1, and for fully incoherent light, M CH approaches infinity and the photo counting noise becoming dominant. The average number of the photoelectrons produced by the detector is computed by: where E½ is the expectation operation, Dt is the sampling period, h is Planck's constant, m is the frequency of the light, and g is the detector quantum efficiency. Background photons are those collected by the sensor which do not originate from the laser transmitter. In most practical scenarios, these photons originate from the sun. This noise can be simulated during LADAR measurement by modelling the number of photoelectrons produced by the background as a Poisson random variable with a mean E½N b given by the following equation [34]: where N b is the number of photoelectrons contributed by the background, including the Poisson noise; S IB is the intensity of background light at the target in units of W=m 2 per lm of electromagnetic bandwidth and D k is the electromagnetic bandwidth in lm of an optical bandpass filter present in the receiver. A B is the target area seen by the receiver. E½N dark is the expected number of electrons contributed by dark current. Figure 12 shows the effect of Photo Counting, Speckle and Background Noise on the received signal. The figure also shows the selected threshold value used to eliminate the effect of this noise. This value is assumed in this simulator to be ten times larger than the maximum noise mean value (which is calculated at a maximum atmospheric transmission, maximum target area seen by the receiver, and minimum target range ð200 mÞ scanned by the LADAR).

Appendix 2: LADAR Range Equation
The range equation is widely used as an analytical tool for computing the power received P r from a target illuminated by a laser pulse containing a given power P t . The LADAR equation is directly analogous to the original RADAR equation [9,23] and can be broken down into several terms that quantify the contribution of the laser propagation elements illustrated previously (transmission, reflection, and reception).
According to [34] the equation is: where P r received signal power at the LADAR detector(watts), P t transmitter pulse power (watts), h t angular divergence of the transmitted beam (rad), Photo-Electrons Threshold Fig. 12 Noisy signal due to photo counting, speckle noise and background noise R ls range between LADAR and the target (meters), q tr target surface reflectance, A tr target surface area (square meters), X tr solid angle of the dispersed radiation (steradian), D r diameter of circular receiver aperture (meters), T a one way attenuation factor, T o optical transmission, 1 intensity reaching the target area from the transmitter (W=m 2 ), 2 ratio of target area to the reflected beam area, 3 intensity at the receiver aperture (W=m 2 ), 4 area for the circular receiver aperture (m 2 ), 5 power captured by the circular receiver aperture (W), 6 attenuation by target reflection, atmosphere and optical transmissions. As mentioned in Sect. 2.3 on page 3, the simulator assumes Lambertian targets. Therefore, X tr , is replaced by the value associated with the standard diffuse targets having solid angle of p steradians [9,23,34].
In the simulation, the laser beam footprint is much smaller than the extent of the target, which means that the target surface intercepts the entire beam (extended targets). This makes the target area A tr equal to the area illuminated by the laser beam, and is given by: With the previous assumptions and substituting Eq. 21 in Eq. 20, the LADAR range equation becomes [3,17,34]: Appendix 3: LADAR Scanning Parameters Table 1 presents the scanning parameters used during performance evaluation for the LADAR simulator.