A Novel Algorithm of Surface Eliminating in Undersurface Optoacoustic Imaging

This paper analyzes the task of optoacoustic imaging of the objects located under the surface covering them. In this paper, we suggest the algorithm of the surface eliminating based on the fact that the intensity of the image as a function of the spatial point should change slowly inside the local objects, and will suffer a discontinuity of the spatial gradients on their boundaries. The algorithm forms the 2-dimensional curves along which the discontinuity of the signal derivatives is detected. Then, the algorithm divides the signal space into the areas along these curves. The signals inside the areas with the maximum level of the signal amplitudes and the maximal gradient absolute values on their edges are put equal to zero. The rest of the signals are used for the image restoration. This method permits to reconstruct the picture of the surface boundaries with a higher contrast than that of the surface detection technique based on the maximums of the received signals. This algorithm does not require any prior knowledge of the signals’ statistics inside and outside the local objects. It may be used for reconstructing any images with the help of the signals representing the integral over the object’s volume. Simulation and real data are also provided to validate the proposed method.


INTRODUCTION
The task of reconstructing the spatial configuration of the sources using their scattered wideband signals received outside the area of the sources location is that of great theoretical and practical interest for various applications. The wellknown tasks of this type include: the optoacoustic detection of inhomogeneities in human tissues (breast tumor detection) [1], and the underground penetrating imaging [2]; a nondestructive analysis of materials [3]. The systems solving these tasks have some common features: (1) the wideband (radar or laser) pulse signal illuminates the object; (2) the scattering object is of a 3-dimensional (3D) shape and composed of point scatters, so the received signal consists of a sum of some scaled and delayed versions of the transmitted signal; (3) the objects which are to be detected are located under a covering surface. The signals from this surface dominate in the dynamic range of the received signals and complicate the process of restoration. Thus, the signals from the surface should be removed. The surfaces in these tasks are the ground surfaces, the surface of the studied material, the skin of some organic body. Among these tasks, the most difficult is the task of medical optoacoustics, since the spatial position of the 3D surface is not known.
Several techniques of "penetrating" imaging are developed in [1,2]. They use different criteria and calculation techniques, based mostly on the idea of cutting off the areas of signals with the maximum magnitude. However, this is not the best criterion. The mathematical technique, using the image gradients' flows for constructing the boundaries contours, has recently become widely used. It consists of building up the contour curve, that satisfies to the minimum of the criterion, in order to adapt it to the boundary of an object. The criteria are various in different works: in [4,5,6,7], the segmentation methods use some special statistical properties of images, which are different in areas divided by the contours. The methods are based on prior knowledge of statistical properties of images and assume a large number of resolution elements in the image. The common features of these most approaches are: the iterative calculating algorithms and the segmentation of the given 2-dimensional (2D) image, when the task of the surface elimination is already resolved or does not exist. The authors of [8] suggested the maximization of the correlation between the ultrasound and MR images for the automatic reconstruction of the 3D ultrasound images.
The paper [9] suggests the algorithm of boundary tracing in the 2D and the 3D images. The boundary is defined as the curve or the surface between the body and the background. The paper [10] develops the program which traces the boundaries of the regions with the definite gray levels in a 2D image, then dissects the boundaries in straight segments end encodes them for compressing the image. The areas restricted by the definite levels of intensity do not necessarily provide the information about the position of the surface, so the algorithms cannot be applied directly to the task of eliminating the covering surface.
Here we address to the optoacoustic task in detail and suggest an algorithm, using the assumption that the objects change smoothly within the inhomogeneities and have the discontinuity of spatial gradients on the boundaries of these inhomogeneities. The algorithm is synthesized to find the lines of the gradients' discontinuities using some mathematical model for these lines. Parameters of this model are estimated by the method of maximum likelihood. The procedure draws 2D (the time index of the received signal, the number of the received signal) curves along which the discontinuity of the signal gradients occurs, removes the areas with the covering surface and leaves the signal areas for the reconstruction of the inhomogeneities. The position of the surface is estimated by a set of the gradients of the signals received along the range coordinate. Then, the detected points of the surfaces are banded in the neighboring signals into the curves, and then, the surface is cut inside these curves. Only then, the restoration of the image is performed. The number of the received signals depends on the characteristics of the receiving aperture and, in practice, may not be very large. Thus, the iterative reconstructing of the active contours may not converge to any reliable result.
The proposed algorithms are investigated by using simulation. The performance of the algorithm is also tested with the help of real signals of the physical model "phantom."

TASK STATEMENT
The task of optoacoustic image reconstruction has the following physical basis [11,12,13,14,15]: the 3D object is placed into some liquid and irradiated by some source. (In our case it is a laser, which generates short pulses, it may also be a radar generating some short high frequency pulses [1].) These irradiating pulses induce an acoustic signal at each point of the 3D object. The acoustic signals from the points are summarized and spread in the 3D space as an acoustic wave. The wave reaches an acoustic receiver, located at some point in the space, and creates some acoustic pressure inside it. This acoustic pressure is transformed into the digital signal in the output of the receiver. If the irradiating laser (or radar) pulse is short enough, the output signal in the receiver has a very high-range resolution. If we have the aperture consisting of a set of such receivers and if the whole aperture covers a large angle of observation, we can restore a 3D image of the irradiated object. If we have a 2D aperture, it gives us opportunity of reconstructing a 3D image. In the case of the 1-dimensional (1D) aperture, looking like a curve, only the integral of the object over the unresolved coordinate can be reconstructed.
Suppose we have N optoacoustic signals Y ( R n , t) (n = 1, . . . , N). According to [11,12,13,14,15], the temporal integral of the acoustic pressure, detected by the transducer, located in point R n , can be described by the following formula: where Y ( R n , t) is the acoustic signal, which is generated by a 3D object when it is irradiated by the inducing source: Here Y ( R n , t) is the integral acoustic pressure in point R n at the moment t, K is the constant proportional to the thermal coefficient of the object volume expansion, exp(−α| r |) is the coefficient of the amplitude attenuation of the signal during its passing through the medium, 1/| r | is the coefficient of the weakening of the wave when it is spread from source O( r ) (the result of resolving the wave equation). O( r ) is in (2) is the shape of the object in the coordinate space r, R n is the vector of the coordinates of the receiver with number n, v is the velocity of the wave spreading (in our case, the velocity of the sound), t is the time index, m( R n , t) is the additive noise in the receiver, which is assumed to be the Gaussian stochastic process, with no correlation between different points R n and the time correlation function ρ n (t) (n = 1, . . . , N), and u(t) is the shape of the laser pulse, inducing the acoustic signal Y ( R n , t). This pulse is very short (∼ 10 nanoseconds in the real system described below). On this supposition, the formula (2) can be simplified as follows (the slowly changing functions can be taken out of the integration sign): If we introduce a new signal X( R n , t) by the formula we will get the following expression for it: Here n( R n , t) is the additive noise with the new time correlation function ρ 1,n (t 1 , t 2 ) (n = 1, . . . , N): If the functions ρ n (t) (n = 1, . . . , N) are narrow enough (i.e., the additive noise in the receiver is closed to the uncorrelated one) we can write a simpler approximation for ρ 1,n (t 1 , t 2 ) (n = 1, . . . , N) as follows: The noise n( R n , t) is uncorrelated between different receivers as before. Exponent α in (3) is generally unknown. The task of its estimation is a separate and a difficult one. In this paper, we will not consider this question, but suppose that α is a priori known. Our task is to get a possibly effective estimate of function O( r ) in the presence of some interfering surface as well as to investigate the quality of this estimating in real conditions.
The function O( r ) is a superposition of the in-question inhomogeneities O obj ( r ) and the surface O sur ( r ), that is, The task of the early medical diagnostics is the detection of small-sized inhomogeneities, that is, the restoration of the image O obj ( r ). The signals from the inhomogeneities have a low amplitude and each of the inhomogenities is located within a narrow time (range) interval. The signal from the surface O sur ( r ) is the signal from the skin and it is generated by a thin irregular curved layer covering a wide spatial range. This signal is very strong and, in fact, it is not zero along the whole time axis (Figures 1 and 2). Each differential element of the surface may not give a significant amplitude of the signal, but a large quantity of such elements, disposed at the identical distance from the receiver, makes a strong contribution into the integral (5). We mean, that the surface spreads into a wide spatial area around inhomogeneities (in a real case, the inhomogeneities can be of several millimeters in a diameter, and the surface-breast skin has an area about a square decimeter). The task of the algorithm is to separate in each signal (5), the areas generated by the surface O sur ( r ) and the object O obj ( r ), and to suppress the areas in signals, generated by the surface O sur ( r ). We will have more convenient conditions for the analysis and the separation of the signals into the areas if we switch to the new coordinate system under a 3D integral (5). Instead of coordinates r x , r y , r z , we will introduce a new coordinate system (τ, ρ 1 , ρ 2 ), where and the coordinates (ρ 1 , ρ 2 ) are disposed in the plane which is orthogonal to the sight line from the chosen receiver. These coordinates supplement (9) to the full 3D coordinates system. Using the coordinates (τ, ρ 1 , ρ 2 ), we can get a new form Now, what we are getting instead of (5) is whereÕ O n (τ) is the new record of the object O( r ) and it presents an integral over the object space in the plane, orthogonal to the sight line from the given receiver. This recordÕ n (τ) is the 1D function of time, and O( r ) is a 3D function. At the same time theÕ n (τ) is an unknown function, different for each new signal X( R n , t), and O( r ) is the function, common for all the signals. Taking (8) into account, we can write We need to find some informative characteristics of the func-tionsÕ n,sur (τ) andÕ n,obj (τ) in (13), which allow to separate the respective signals. We can suggest the time derivatives of these functions as the informative characteristics. These derivatives have their maximums (of absolute values) at the boundaries of the object (at the front edges ofÕ n,sur (τ), O n,obj (τ), and at the back edges of these functions, resp.). At the edges, these derivatives are close to delta functions. Anyhow, this is true about the inhomogeneities with the shape close to the spherical one (with a small radius) and for the surfaces of some arbitrary shape and size, but thin, however. Very often, the task of the medical diagnostics has the similarity to the task of detecting a small-sized inhomogeneity of a spherical shape. We consider the time-derivatives of the signals given by (14) and design them as Gr( R n , t). Using (14), we can write wherem n (t) is the additive noise with the new-time correlation function. This correlation function can be calculated directly and it equals to ρ 2,n (t 1 , t 2 ) = ∂ 2 ρ 1,n (t 1 , t 2 )/∂t 1 ∂t 2 (n = 1, . . . , N). All the noisesm n (t) are uncorrelated between the different receivers, because the transformation (15) is being performed independently between the different positions.
We can easily see that (15) can be replaced by Now we can formalize the problem of signal separation. Further, we will search for the function dÕ n (τ)/dτ as a sum of a certain slow function and an unknown number of delta functions with some arbitrary amplitudes and location of maximums Here A 0n (τ) is the slow function and δ(τ) is the delta function.
Parameters I n , A in , and τ in and the function A 0n (t) are unknown and should be estimated. The approximation (17) assumes that the form of the signalÕ n (τ) along the range τ is a smooth function of τ except for some areas, where the inhomogeneities and surfaces are located; and the derivatives dÕ n (τ)/dτ have the discontinuities on the edges of these areas.
This approximation does not fully correspond with the physical properties of the signals, of course. But, the approximation (17) permits to extract the delta-form peaks in the derivatives of signals and to detect the local objects with using asymptotic methods [16]. A method of estimating parameters I n , A in , and τ in , and the functions A 0n (t) is given below in Appendix A.

FULL ALGORITHM OF IMAGE RESTORATION UNDER THE SURFACE
Formulas (A.11) and (A.13) give the estimates of parameterŝ A in ,τ in , andÎ n (i = 1, . . . ,Î n ; n = 1, . . . , N); overall, the algorithm of building and using the separating curves consists of the following operations.
(2) The construction of the curves of the gradients' discontinuity. The curve with the number i = i 0 is a set of parametersτ i0n for a certain number i = i 0 and for all the numbers n (n = 1, . . . , N), constructed on the basis of the whole set of the received signals. This curve T i0 = (τ i01 ,τ i02 , . . . ,τ i0N ) can be considered the boundary of the local object and, thus, it can be used as the line separating the signals into the areas. If, in addition, this region is characterized by the maximum values of the estimates |Â i0n |, it can be considered exactly the area where the signals from the surface are located.
the values of the signals within this curve should be set to zero. If the surface lies between the receives and the unclosed curve T i0 = (τ i01 ,τ i02 , . . . ,τ i0N ), then we have to set all the signals at axis t in the intervals (0,τ i0n ) (n = 1, . . . , N) equal zero. If the surface lies behind the inhomogeneities along the range, then we have to set all the signals at axis t at the intervals (τ i0n , T) (n = 1, . . . , N) equal to zero (here T is the last time point of all the received signals). (4) After this operation, we can apply the image reconstruction procedure described in [17]. This procedure comprises two operations (in a case of the 2D restoration).
(a) The summation of all the signals in the plane of the image reconstruction performing the transition from the time coordinates to the spatial coordinates of the image: (b) We will design the 2D Fourier transform of (18) as F Z ( ω), where ω is the variable of the spatial frequencies. (c) The multiplication of F Z ( ω) by the filtering function H( ω): where τ pulse is the length of the inducing pulse u(t).
It should be noted that formula (19) was exactly derived in [17] only for the Gaussian form of the pulse The filter (19) suppresses the low frequencies down to zero, retains the middle frequencies without any changes, and suppresses the high frequencies; (d) The reverse Fourier transform of the result received by multiplying gives the final estimation of O obj ( r ). It is clear from formulas (19) and (20), that the essential parameters of the algorithm are the velocity of the wave propagation v and the length of the inducing pulse τ pulse .
It should be noted that there are two options for the implementation of the algorithm in constructing the curves T i0 = (τ i01 ,τ i02 , . . . ,τ i0N ).
(A) By the analytical calculation of (A.10) and its maximization. (B) By using the interactive computer work mode. In this case, we have to take into account the following considerations: X( R n , t) is a function in 3D space; R n is the point of the aperture where exactly the receivers are located (e.g., a semisphere [1] or a plane [18]), t is the time axis for the signal. We can assume that the receivers are located in a single-plane layer, for example, along a certain curve in the plane XY. This assumption retains the applicability of the technique for any 3D shape of the aperture, since for each new layer (along the Z-axis), we can use the procedure again. In case we have a rather large receiving aperture with the receivers located closely to each other, the signals (15) and (16) will vary continuously between the receivers. Thus, the processing should include the following operations: (1) to reconstruct on the display of the computer all the N modules of the signal gradients: received within the single plane (Figures 2 and 3); (2) to set to zero all the signals on the left-hand side or on the right-hand side (depending on the specific location of the surface) of the curves T i providing the maximums to the values of (21). In the interactive mode, the positions of these curves should be indicated by an analyst with using the "mouse." Below, we will discuss this technique and demonstrate the procedure.

TESTING THE ALGORITHM BY USING SIMULATION
The computer simulation model of the signals is useful for testing the performance of the algorithm. All the objects (the four spheres of different diameters and the interfering surface) were simulated by using "OpenGL" package of 3D graphics [19]. The surface model is a set of polygons simulating a certain large sphere. All the polygons are equally thin (about a diameter of the smallest sphere).
The number of the receivers is 32. They are arranged along the circle with a radius of 60 mm in plane XY and cover the observation angle of 120 degrees. Figure 4 shows the whole true object in plane XY, where the receivers are located and it is the area of the image to be restored as well.
Each position receives a signal at the time interval of 134.228 nanoseconds. The number of the points in the signal is 596. The velocity of the sound is 1500 m/s. The signal covers the range interval of 120 mm. This interval was taken as the size of the volume under investigation. The arrangement of the receivers is shown in Figures 5, 6, 7, and 8.
The signal (14), received by the position under number 17 (in the center of the receiving aperture) prior to cutting, is shown at Figure 1 as the function of the range. Figure 2 presents the set of the magnitudes of the gradients of 32 signals, calculated with the help of formula (21). The signals (14), which are the signals received from the four spheres and the surface were also simulated and computed in the "OpenGL" package. In Figure 2, the (ρ = tv)-axis of ranges is horizontal and the n-axis is vertical. The area (to the left) occupied by the surface is rather distinct. The surface is exactly between the receivers and the spheres and it simulates the breast skin. This is the area of the maximum values of the signals (14) and the maximum values of the signal gradient magnitudes (21). In general, the surface covers almost the whole plane (ρ = tv, n), but in the middle and on the right-hand side area in Figure 2 the levels of the signals and the gradients from the surface are much lower. That is why, we may cut off only the maximum values on the lefthand side area of Figure 2. The cutting line was drawn by the mouse in the interactive mode and recorded at the operative memory. After that, all the 32 signals on the left-hand side area of the curve were set to zero. The result of the cutting operation is shown in Figure 9. We can see from Figure 9 that the signals from the spheres and from the part of the surface overlapping with the useful signal are retained in the plane (n, t = ρ/v). Figures 5, 6, 7, and 8 present the reconstructed images of the four spheres: Figures 6, 7, and 8, under the surface and Figure 5, with no surface at all.
As it was said, all the signals cover the range interval of 120 mm (beginning from the range which is equal to zero). So the volume within which the restoration of the image is principally possible, has the dimensions of 120 × 120 × 120 mm. As our aperture has only 32 receivers, located in the plane, the 2D space of the image restoration is 120×120 mm, that in pixels equals to 596×596. The central point of the im- age frame has the range of 60 mm from the central receiver. The scale (in mm) is shown along the axes X and Y in all the pictures. The arrangement of the receivers is shown at the bottom of the figures. Figures 5, 6, 7, and 8 present the result of the image restoration using the algorithm [17]. The image is shown in the plane XY. Figure 5 is the restored image of the four spheres without any interfering surface (only spheres). Figure 6 presents the result of the image recovery with the surface present, when only the summing up procedure of all the signals is performed in the plane of the image (the first stage of the algorithm [17]). Figure 7 shows the result of the image restoration under the surface after the optimal space filtration (the second stage of the algorithm [17]). Figure 8 demonstrates the restored image after the process of the surface cutting algorithm and procedures of summing and filtration. The level of the surface has become lower,  Figure 9: Magnitude of the gradients of all the signals after cutting off the surface. and the resolution of each of the spheres is improved. The smallest sphere placed at the greatest distance from the receivers can be observed almost as sharply as in Figure 5 (only spheres).
All modeling was performed without taking into account the noises in the receiver. To evaluate a comparative efficiency of the described algorithms, some calculations of the potentially reachable signal/noise ratios are given in Appendix B.

TESTING THE ALGORITHM BY USING THE REAL SIGNAL FROM THE PHANTOM
The real optoacoustic system with the arc-array transducers processing the optoacoustic signals was described in detail in [20]. The aperture has 32 rectangular receivers of 1.0 × 12.5 mm dimensions, and the distance of 3.85 mm between them. The transducers are located on the circle with the radius of 60 mm. The real physical model was a sphere with the diameter of 0.8 mm, placed in milk. The milk was diluted with water to obtain optical properties of the medium close to the ones of the breast tissue. The optical absorption coefficient of the sphere was about 1.0 cm −1 . This value is typical of some light absorption in tumors [20]. The sphere is disposed in the near zone, approximately above the central receiver, at the distance of 19 mm from it.
The laser radiation comes along the Y axis. The energy of the laser pulse is within the range of 0.025-0.050 J to comply with the regulations for the medical procedures, which require that the density of laser radiation at the surface of the breast should not exceed 0.1 J/cm 2 . All the receivers are arranged equally and they cover the angle of 120 degrees. Each position receives the signal with the rate of 66.667 nanoseconds. The number of the points in the signal is 1200. The range interval covered by the signal is that of 120 mm. This interval was taken as the size of the volume to be investigated. The arrangement of the receivers is shown in Figures  11 and 12. Figure 3 presents the set of the gradients' magnitudes of all the 32 signals, calculated by using formula (21), for all the real signals. The strongest part of the surface has been already cut off from the signals previously and, thus, is not shown in Figure 3. However, the significant elements with the surface areas still remain. We can see in Figure 3 that there are several areas (on the right-hand side and in the middle of the picture) occupied by the surface. These are the areas with the maximum values of the signal gradients. Several lines and several areas for signal cutting are distinctly visible. The brightest area on the right-hand side and in the center of the picture was the first to be cut off in the interactive mode. Then, on the left-hand side of the picture, a new bright area stood out, that was cut off as well. The final result of cutting is shown in Figure 10. We can see that only signals coming from the sphere and some background noise remained in the plane (n, t = ρ/v).
In Figure 11, we present the image, constructed in the plane (X, Y ), where the receivers are located, and prior to cutting off the surface-related signals. The image was reconstructed in the frame of 120× 120 mm or 1200 × 1200 points. The recovered image is the result of the summing and filtration, performed according to [17]. Figure 12 shows the restored image after removing the surface. We can see that, in fact, the sphere only remained in the image.

DISCUSSION
The proposed algorithm makes it possible to reconstruct the edges of local objects and the boundaries of the surface covering these objects. The data used in the algorithm, are the spatial gradients of the received signals. This method permits to reconstruct the picture of the surface boundaries with a higher contrast than that of the surface-detection technique based on the maximums of the received signals. This algorithm has also an advantage over the method of the active contour; it does not require any prior knowledge of the signals' statistics inside and outside the local objects, and it does not function as an iterative procedure either. This algorithm may be used for reconstructing any images with the help of the signals representing the integral over the volume of the object (5), but as for the optoacoustic signals, it has already been tested on the digital model and real signals. Figures 2 and 3 illustrate that the signal gradients' magnitudes (21) are good indicators for localizing the surface and detecting the inhomogeneities in the volume. The procedure using the complete set of signals for determining the area occupied by the surface is suggested. The algorithm constructs the curves t(n) showing discontinuities of the signal derivatives (the time index of the discontinuity is t = ρ/v, where ρ is the range value in the figures, the number of the signal is n). These curves t(n) can be drawn by using the mouse in the interactive mode. Figures 8 and 12 illustrate that the process of cutting off the area occupied by the surface, this leads to improving the images of small inhomogeneities.  Figure 11: The recovered image of real phantom (one sphere in milk medium) prior to cutting off the signals from the surface. Figure 12 shows that the algorithm can be applied to the real experimental system [20] with the real signal energy and the real contrast levels.
It is useful to discuss the computational demands on the proposed algorithm.
The main computational requirements are imposed to the procedure of the image restoration, that is, to the operations, described by (18), (19), and (20). The operation of the signals summing (18) in the window of 1200 × 1200 pixels was calculated within 52 seconds at the computer with 256 Mb RAM and 1300 MHz clock rate. The operation of filtration (19) took 16 seconds. The procedures of the surface elimination are less laborious; searching for local maximums and bunching them into the curves took 7 seconds with using 32 signals each of 1200 points of the length (formula (A.10)) with the approximate calculation of integrals along the time. Performing this procedure in the interactive mode is slower, but it is more reliable. (17) We chose the approximation (17) to insert it into (16) to locate the peaks in the gradients (16). The parameters of these peaks, τ in , A in , and their number I n are unknown and must be estimated. In (17), τ in are the time indexes of the local edges ofÕ n (τ), and the sign of A in is the sign of the gradient at the local edge ofÕ n (τ). I n is the number of discontinuities (proportional to the number of the separate local objects).

A. ESTIMATING ALL THE PARAMETERS OF
It should be noted that the delta function δ(τ) is a generalized function with the property of filtrating the single point of the function under an integral, that is [21], By inserting (17) into (16) and by using the property (A.1), we will get An approximation (17) leading to formula (A.2) is the mathematical assumption and of course does not always correspond with the real physical conditions. The proposed model (17) is only an asymptotic approximation to the real physical model. But the digital modeling and processing of the real signals show (in the sections describing the testing algorithm) acceptability of such approximation. In this section, we will estimate the parameters of this model (A.2) by the method of maximum likelihood.
As it was said above, we assume that the pulse u(t) is very short compared to the interval of constancy of the slow function A 0n (t). On this assumption, it is easy to see from the (A.2) that the peaks of all the gradients in the signals have the width as the width of the inducing pulse u(t). In other words, the edges of the inhomogeneities can be detected with the accuracy not exceeding the range resolution of the given system.
We have N signals (A.2), and our task is to make the estimations of all the unknown parameters A 0n (t), A in , τ in , and I n . The most important parameters are I n and τ in (i = 1, . . . , I n ; n = 1, . . . , N). These parameters describe the shape of the curves which separate the local objects. The curves allow to detect and remove the strongest signals from the surface and to find all the other local objects of smaller sizes.
Further, we will assume that the temporary correlation function of every signal (1) ρ n (t) (n = 1, . . . , N) is narrow enough, so the noises in all measurements of the signal can be considered statistically uncorrelated. In this case, the second derivative of the function ρ n (t) (n = 1, . . . , N) will also be a narrow one and approximately of the same duration as the function ρ n (t) itself, and all the measurements of gradients (15) and (16), having the time correlation function ρ 2,n (t 1 , t 2 ) = ∂ 2 ρ 1,n (t 1 , t 2 )/∂t 1 ∂t 2 (n = 1, . . . , N) (see (6), (7)), can also be approximately considered uncorrelated in time. On this assumption we can write the logarithm of likelihood function LnP [22] for the functions Gr( R n , t) (A.2) under the specific values for parameters A 0n (t), A in , τ in , and I n as follows: Here, LnP n is a logarithm of the likelihood function for signal gradients in the pulse with number n: In (A.5), a designation was introduced: It can be seen from (A.4) and (A.5) that maximization of the whole LnP breaks up into the independent maximization of each of the functions LnP n . The maximization of (A.5) over the parameter A in can be performed exactly. This maximization gives the following estimations:Â Here C 0 is the energy of the pulse: The insertion of (A.7) into (A.5) gives the new view of the function LnP n depending on τ in and I n only as follows: The first term of (A.9) does not depend on τ in and I n . So, we have a function LnP (1) n for maximization on these parameters: (A.10) The likelihood function (A.10) is analogous to the likelihood function in a process of detecting the radar targets in a radar receiver having a square detector when the number of targets I n is unknown [23,24,25]. In our task, the edges of local objects perform a role of targets in the space of gradients. In a correspondence with [23,24,25], this detection and I n estimation should be performed by the following algorithm: with no prior knowledge of the target number I n and their position τ in , we have to use a maximally possible range (0, T) of values τ in (defined by the experimental conditions) and to construct the likelihood ratio: 11) in every point τ in of the whole range.
Formula (A.11) is the likelihood function for the local edge in the point τ in . It can be seen from (A.11) that the local likelihood equals the ratio signal/noise for the parameterÂ in , where σ 2 noise is the dispersion of the noise in the signal gradient function. It can be expressed through the energy of the pulse as follows: (A.12) After the ratio (A.11) is formed, we have to check a condition of exceedingÂ in over the noise, that is, we have to check the next condition in every point τ in : We make a decision about the new detected maximum in the signal gradients if (A.11) exceeds some threshold. In statistical measuring tasks, a value of the threshold is often taken in an interval 1-9. The total number of maximumsÂ in , satisfying (A.13), gives the estimateÎ n of all the front and back edges in all the local objects detected in the signal under number n.
The time positions of these maximums are given by the values τ in in (A.11). Now, it should be mentioned that (A.10) and (A.11) comprise the unknown functions A 0n (t) inside the function S n (t) (formula (A.6)). It is natural to assume that |A 0n (t)| |A in | (i = 1, . . . , I n ), that is, the boundaries have the higher contrast and they are more visible in the space of the gradients than the smooth parts of the derivatives. In this case, we can put A 0n (t) = 0 in (A.10) and (A.11) (as the first step of the calculations at any rate), and the likelihood function for maximization LnP (1) n will obtain the following form: An algorithm of getting the estimates of the slow background A 0n (t) is described below. After these estimatesÂ 0n (t) are obtained, we have to use the function for LnP (1) n in the view (A.10), but formula (A.14) may be used as the first approximation. The sense of the maximization of (A.10) or (A.14) is obvious; the best estimates of τ in and I n provide the maximum for the correlation of the gradients (16) (after leaving the slow background A 0n (t) out of the gradient Gr( R n , t)) with pulse u(t). If u(t) is a short pulse, the maximization of (A.10) or (A.14) simply leads to the search of all the maximums of Gr 2 ( R n , t) along the time axis.
The last step is the evaluation of the slow component of the gradients, that is, the functions A 0n (t) (n = 1, . . . , N). When all the valuesτ in andÎ n are obtained (i = 1, . . . ,Î n ; n = 1, . . . , N), the expression for LnP n will have view (A.10) with inserted estimatesτ in andÎ n into it.
By maximizing (A.10) regarding S n (t), we will obtain the equation for S n (t) as follows: The first approximation for the solution of (A.15) is located in the vicinity of A 0n (t) = 0, and for the short pulses, u(t) has a view: Gr R n ,τ in u t −τ in , (A. 16) where Strictly speaking, we should return to the operation (A.11) after getting (A. 16) and repeat all the calculations again including (A.15). In other words, the process of the simultaneous estimation of the background A 0n (t) and the parameters τ in and A in must be iterative. In a case of the low levels of A 0n (t), the single iteration will be enough. Now it is necessary to say that the described algorithm gives estimatesτ in with accuracy equal to a discrete Del of the data receiving signals. The more accurate measuring of the gradients maximums position will demand the more accurate evaluation τ in in the functional (A.10). It is possible to use the accurate methods analogous to the radar methods of the target location measurement. But the image can be restored only with the resolution Del, even in the absence of the interfering surface. So, the determination of the edges of local objects with a higher accuracy is not necessary, but may appear more labor intensive.

B. COMPARISON OF TWO ALGORITHMS IN SNR
The signal from a sphere S sph (r) as a function of the distance from the receiver r can be described by formula [13,17] as follows: S sph (r) = π Rad 2 sph − r − R 0 sph 2 , (B .1) if Rad 2 sph ≥(r − R 0 sph ) 2 and S sph (r) = 0, otherwise, where R 0 sph is the position of the sphere's center, Rad sph is the radius of the sphere. Further, we will suppose that the surface is a sphere, which is empty inside and has the thickness of its sheath equal to ∆ sur .
The signal from the surface S sur (r) can be described by formula [13,17] as follows: if Rad 2 sur ≥(r − R 0 sur ) 2 and by S sur (r) = 0, otherwise. Here, R 0 sur is the position of the surface's center and Rad sur is radius of the surface sphere.
The full signal S(r) in the receiver got from the distance r will be equal to S(r) = S sph (r) + S sur (r) + n(r), where n(r) is the Gaussian noise uncorrelated between the neighbor data measurements with a dispersion σ 2 noise . A ratio of the sphere signal to the sum of the noise and the surface signals (SNR) in (B.3) is equal to Q sph/(sur +n) (r) = S 2 sph (r) S 2 sur (r) + σ 2 noise .

(B.4)
Now, we will compare this SNR got after the surface eliminating by two methods: (1) cutting off the maximum surface signal; (2) cutting the surface in the gradients space.
(1) If we cut off the maximal level of signal up to the level of the first neighbor minimum, we will receive the new maximal signal of the surface S sur,1 min approximately as follows: