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:


Tools for Teaching Regression Concepts Using Dynamic Graphics

This paper extends work on the construction of instructional modules that use graphical and simulation techniques for teaching statistical concepts (Marasinghe, et al. 1996; Iversen and Marasinghe 2001). These modules consist of two components: a software part and a lesson part. A computer program written in LISP-STAT with a highly interactive user interface that the instructor and the students can use for exploring various ideas and concepts comprises the software part. The lesson part is a prototype document providing guidance for instructors for creating their own lessons using the software module. This includes a description of concepts to be covered, instructions on how to use the module and some exercises. The regression modules described here are designed to illustrate various concepts associated with regression model fitting such as the use of residuals and other case diagnostics to check for model adequacy, the assessment of the effects of transforming the response variable on the regression fit using well-known diagnostic plots and the use of statistics to measure effects of collinearity on model selection.

Instructional Modules
A primary goal of a new NSF-supported project targeted at improving the effectiveness of undergraduate statistics courses has been to develop instructional tools that are easily adaptable for general use independent of specific courses. These tools will provide statistics instructors with the capability of supplementing their teaching methods by presenting and illustrating statistical concepts more effectively than is possible using conventional instructional methods. They allow exploration of important concepts using graphical displays that are easy to use and provide instantaneous visual feedback, thus encouraging active learning.
Under this project, we have developed a number of instructional modules designed to illustrate various statistical concepts and provide important insights into the application of these concepts. Use of these modules in teaching will enable students to get more meaningful learning experiences than are otherwise possible from traditional instructional methods alone. Exercises designed to accompany these modules will further reinforce these learning experiences. The software component of these modules employs a combination of state-of-the-art computing hardware, a statistical programming environment, high resolution color graphics, computer simulation, and a highly interactive user interface. The modules can be used both as classroom demonstration tools and as self-paced exercises for individual students or small groups. The instructional component is a prototype lesson that includes a description of concepts to be covered, instructions on how to use the module, and some exercises with solutions. Course coordinators and instructors can modify or augment these lesson plans to create their own class assignments.

Implementation of the Software Components
The software components of the modules described here have been programmed in the Lisp-Stat language. The system provides a computing environment with the following important features: It is a powerful object-oriented programming language that allows rapid prototyping and development.
It allows development of modules using true dynamic graphics.
It allows links to routines written in Fortran or C.
It can be installed on the different instructional computing platforms in use at many educational institutions (Unix workstations, Macintoshes, and Desktop PC's running Windows software).
It has a large core of statistical functions available.
Another important advantage of using Lisp-Stat for developing software is that, for noncommercial applications, it is available without cost. For a good introduction to Lisp-Stat and other Lisp-Stat based software, see de Leeuw (1994).
We have developed user interfaces that are consistent across modules. Students need only execute one command to start a module. All further interaction is through a mouse-driven interface. One may click on a push button to initiate action, click and drag on a menu button to select items from a pull-down menu or click and move slide-bars to set values of various parameters. Some of the dynamic graphical techniques used in the modules are discussed in detail in Cook and Weisberg (1994).

Lesson Prototypes for Instructional Modules
For each instructional module, we have developed, in LaTeX format, a sample instructional document. These documents consist of An introduction to and a short description of the important statistical concept(s) to be covered.
Objectives for the instructional module.
Instructions on how to execute the software component of the module.
Warm-up exercises to familiarize the student with the software module.
Formal exercises and questions requiring execution of the software(usually following a reasonably precise set of instructions), careful thought, interpretation of results, and explanation of conclusions.
Notes for the instructor, including comments on what is expected to be understood from doing the formal exercises.
Instructors can either use these documents or customize them for their own purposes and classes. The documents (perhaps without the notes for the instructor) could be made available to the student as a classroom handout or online. Becker, Cleveland, and Wilks (1987) review many of the dynamic graphic techniques such as linking and selecting. Many of these techniques have been extensively used in the exploration and the analysis of multivariate data (see Cleveland and McGill (1988)). Uses for rotation and animation techniques in regression diagnostics were first developed by Cook and Weisberg (1989), and later a more elaborate presentation appeared in Cook and Weisberg (1994).

Related Work
However, until recently, little use of these techniques has been made for instructional purposes though the possibilities are many. Saunders (1986) describes the use of dynamic graphics to produce "moving visualizations designed to introduce difficult concepts, reinforce mathematical ideas, and explore the techniques of probability modeling." The visualizations were used in a distance-learning television program produced by the BBC. The paper, for example, describes visualizations to illustrate the effect that parameter changes have on binomial, Poisson and bivariate normal distributions.
Another example is the STEPS project, a UK consortium based in Glasgow involving nine departments in seven universities, which was conceived for the purpose of developing problem-based teaching and learning materials for statistics. Modules were developed around specific problems in several subject areas and incorporated computational and graphical tools to assist in the exploration of the statistical ideas encountered in solving these problems. The problem-based approach was used to motivate interest in the students by selecting problems in their own areas of study and also because it allowed integration with other more standard laboratory materials. Trumbo (1994) uses graphics and simulation to illustrate elementary probability concepts, using programs written in QuickBasic and giving rough equivalents in Minitab.
Even fewer examples appear in the area of regression where the techniques and concepts are ideally suited for illustration using dynamic graphics. Apart from the introductory examples in Tierney (1991) and the work of Cook and Weisberg (1994), a paper that describes software constructed for this purpose is Nurhonen and Nurhonen and Puntanen (1992). Anderson and Dayton (1995) present program code written to demonstrate various features available in the Lisp-Stat language. Rather than providing a complete set of educational modules, they illustrate how this language can be adapted for building instructional tools to enhance the teaching of various concepts in regression. That approach may not be suitable for instructors averse to learning the Lisp-Stat language at the required level for developing modules. On the other hand, no programming ability is required to use the modules presented in this paper. In addition, the lesson plans provided serve as templates for instructors for creating their own lessons. Users of these modules may find different ways to use them other than those suggested in the paper or the lessons. The present modules are also extensible in that anyone may add a module to the system without affecting current users of the modules.
A commercially available multimedia package ActivStats written by Paul Velleman contains several interactive modules that illustrate various concepts in regression. However, these modules are built around an introductory statistics course, and instructors may find integrating ActivStats with a regression course difficult. Also, ActivStats is not free.
There are several archives of JAVA applets available for regression. An example is the VESTAC system at www.kuleuven.ac.be/ucs/java/, described in Darius, Michiels, Raeymaekers, Ottoy, and Thas (2002). While some of these JAVA applets may be useful for demonstrating the same regression concepts covered in this paper, no lessons are provided with the applets. Using the JAVA applets requires a network connection to access the applets. This can result in slower execution and a less responsive interface than locally installed modules. Modules described in this paper are based on Lisp-Stat and can be customized using the Lisp-Stat language; however, the JAVA applets cannot be customized.
A more recent addition to the literature is the Cook and Weisberg (1999) regression text accompanied by the software package Arc, also written in Lisp-Stat. Although designed specifically for performing analyses described in the book, Arc could conceivably be used independently to demonstrate selected regression concepts. A more useful suggestion for instructors using the above text and Arc, is to use the modules described in this paper as a supplement to their course.

The Regression Modules
In this article we describe a set of instructional modules that have been specifically designed to aid instructors in teaching introductory and advanced regression methodology. For simplicity of use and ease of construction, the software component consists of modules each covering a different but interrelated set of regression concepts. The five regression modules are described in detail in Sections 2.2 through 2.6. Section 3 describes how to obtain the software described in this section. Readers may find it helpful to install and use the software while following the descriptions in Sections 2.2 through 2.6.

Regression Concepts
Statistics students from different disciplines often have various levels of ability in mastering some of the more sophisticated modern regression techniques they learn from textbooks and lectures. Interpretation of results produced by standard software packages such as case deletion statistics or residual plots is a complex task and can lead to confusion and improper use of such statistics. An adequate understanding of the concepts behind these techniques and some experience in using them can help alleviate this problem. However, students enrolled in regression courses cannot acquire the necessary experience through involvement in real statistical data analysis projects alone because these projects are often too time consuming to incorporate more than a few into a course. Some examples of the kinds of concepts that students need to understand for developing an ability to interpret regression computations are: Different graphical displays highlight different relationships among variables. To explain the relationship between a Y and an X variable in a multiple regression model, both Y and X must be adjusted for effects of other explanatory variables in the model.
Interpretation of case statistics from a regression analysis is often far from straightforward. Not only does the presence of more than one extreme observation tend to complicate their interpretation, different case statistics provide fundamentally different types of diagnostic information.
If diagnostic plots indicate departures from assumptions (e.g. nonnormality of residuals, nonhomogeneous errors, etc.), a transformation of either the Y or the X variable (or both) may result in variables that more closely match the model assumptions.
While standard homework assignments, lab exercises, and class projects provide mechanical practice, illustrate a few ideas, and help to introduce these concepts to students, true understanding comes only after extensive experience. Repeated interaction with graphical displays of simulation results allows students to be exposed to pseudo-experiences that will facilitate their understanding of important concepts. For example, by plotting straight lines fitted to regression data generated from a known fixed model repeatedly, the student can understand the difference between fitted models and the true model. In another example, students can dynamically change the value of a data point to study the effect of that point on the fitted regression line. By controlling the way that a change is effected (e.g., by holding X fixed and changing Y in simple linear regression), the program can highlight a single property of a case statistic (such as the simple fact that leverages depend only on the predictors). Examples like these illustrate the potential benefit of incorporating modules like ours into a traditional regression course.

Exploring Least Squares Fitting
The regteach1 module is useful for illustrating some of the fundamental concepts related to simple linear regression. For example: A single summary statistic like a correlation coefficient or R 2 , by itself, cannot be used to interpret the strength of a relationship. A scatterplot is an essential component of examining the relationship between two variables.
It is important to understand the idea of least squares fitting. It can be demonstrated that one may not always be minimizing the sum of squared deviations when "fitting a line by eye".
Magnitudes of the residuals from a regression depend on the fitted line. Thus a simple residual plot can reveal a lot about the goodness of the fit.
The frame on the left-hand side of Figure 1 shows the initial view of the regteach1 module window. Changing the slope and intercept values in the slide-bars will dynamically change the slope and intercept of the plotted line and update the numerical coefficients in the box above the plot. Pushing the Select Data menu button and using the resulting pull-down menu allows the user to select a data set from a list of simple regression data sets. The frame on the right-hand side of Figure 1 shows the regteach1 module after the OAK-SEEDLING data set has been selected. Pushing the Residuals button will produce an additional window containing a plot of residuals from the current fitted line. The signed deviations are displayed as line segments drawn from a zero baseline, plotted at the corresponding x-values. The residuals will change as the slope and intercept are changed in the main window. This can be used to demonstrate dynamically the dependence of the magnitudes of the residuals on the fitted line (such as the fact that as the fit improves the residuals get smaller), as well as identifying various patterns in the residuals (such as curvature) for diagnosing the fit. Figure  2 shows the residual plots corresponding to two different lines fitted to a data set. The user can attempt to find the best fitting line to the selected data "by eye" (graphically) using the sliders to change the slope and intercept; the plotted line will be updated dynamically. When satisfied with the line fitted "by eye", the user can push the Fit Least Squares button to display the least squares line fitted exactly to the data and the resulting summary statistics (including regression coefficients, sum of squared residuals, and the residual standard deviation) as shown in Figure  3. Comparing the least squares regression line to several different "by eye" fits provides the user the opportunity to learn the meaning of the least squares criterion for fitting a line to data. The Fit Least Squares button also creates a scatterplot of the residuals plotted against the X variable for making a visual assessment of the fit (e.g., curvature exhibited by the residuals may indicate nonlinearity in the data). The Correlation button will display the between X and Y. In addition to a graphical fit, a user may also attempt to fit the best line numerically, by attempting to minimize the value of the RSS, or better, by using the RSS "thermometer" displayed on the right margin of the primary window, both activated by pressing the Residual SS button.

Properties of Regression Coefficient Estimates
The regteach2 module is constructed to illustrate some of the concepts related to statistical properties of parameter estimates in the simple linear regression model. In particular, students will have an opportunity to examine the sampling distribution of the slope estimate dynamically. Some of the ideas that can be explored with this module are: The fitted line passes through the "center" of the data, i.e., through the point .
Variability in the data affects the accuracy of estimation of the regression coefficients.
The estimate of the regression slope is symmetrically distributed around the true value of the slope.
Parameter estimates depend on the data only through a few summary statistics.
If it is assumed that the errors are a random sample from N(0, ), the slope estimate b is distributed normally with mean , the true value of the slope, and variance . Thus it is easy to see that for fixed values of the regressor variables, the distribution of the estimate depends only on .
The initial window of the regteach2 module is shown in the top frame of Figure 4. A value for the variance is first selected for generating data from the model (displayed in the box near the top margin of the plot) by moving the slide-bar. For each push of the Simulate Lines button, seven new y values corresponding to the seven x values specified are generated, and the straight line fitted to the data is plotted. By repeating this procedure one can observe dynamically the variation in the lines fitted to different data sets generated from the same model. The slopes of fitted regression lines constitute a draw from the sampling distribution of the slope estimate b. By pushing Dist. of b prior to starting the simulation process, the user can observe these values being accumulated in a dynamically updated histogram. As the histogram is updated with more and more samples of the estimated slope, the user is able to visualize that the distribution begins to take the appearance of a Normal distribution centered at 2, the true value of the slope. The Save 15 Lines button, when pushed, will plot 15 such lines on the same graph for a fixed value of and save a copy of the graph in a separate, smaller window. This is useful for comparing sets of lines obtained for different choices of , as exhibited in Figure 5. These graphs clearly demonstrate that the variation in the slope of the lines fitted to data is dependent on the actual variability in the data.

Diagnostic Use of Case Statistics
Many regression programs in standard statistical packages produce a large number of the regression diagnostics proposed in the literature. However, not much attention is paid to describing or indicating how such statistics should be used or interpreted. Rather than providing students with additional tools for judging the adequacy of a fitted model, this vast array of statistics has added to the confusion of many students.
The regteach3 module has been designed to illustrate important concepts related to the use of some of the more popular case statistics in regression analysis: The leverages depend only on the explanatory variables. Cases away from the centroid of the data have large leverages compared to those near the middle.
The predicted response for cases with high leverages is largely determined by the observed response. Hence the fitted line is constrained to pass as close to the corresponding data point as possible.
Deletion of a single point can have a large effect on the fit of a model. A case is determined to be influential if its deletion substantially affects the parameter estimates.
Cook's D is a combination of an outlier measure and a leverage measure. An influential observation (i.e., one with a large D value) that has small leverage is an "important outlier." When regteach3 is started-up, the main window (top left frame in Figure 6) shows a scatterplot of a working data set along with an overlaid regression line (displayed in magenta; the regression statistics corresponding to the fitted line are also displayed in the same color). The three secondary windows display, clockwise from the main window, indexed plots of the studentized residuals, Cook's D statistics, and the leverages, respectively, computed with respect to the fitted regression line. Initially, all four windows of the module are in selecting mode; i.e., a point may be selected by clicking the mouse pointer near it. By pushing the Point Moving Mode button, the main window can be changed to the point-moving mode. In addition, all four plots are linked with each other, so that if a point is selected in one of the secondary windows, the serial index corresponding to the point will identify the point in all four windows. In the point-moving mode, points in the scatterplot can be dragged into other positions; however, in regteach3 they are constrained to move only in the vertical direction, so that the value of the X coordinate of the point moved remains fixed. To demonstrate a case statistic like leverage, we will first move to a secondary window and select one of the points indexed, say, 18 or 19. The selected point will be highlighted and labelled. These points are furthest from the "centroid" of the x values (the centroid here is the mean , correspond to the two largest leverages, and lie above the 2p/n reference line shown on the leverage plot. Now, moving to the main window, the selected point is dragged to a new location (i.e., the point will have a new Y coordinate while the corresponding x coordinate remains the same). The regression line will be recomputed and plotted in a new position. The original line will still appear on the graph (in a light green color) for comparing with the refitted line. Notice that this action leaves the leverage plot unchanged, while the other two plots, as well as the estimates of the parameters, are all updated dynamically to correspond to the line fitted to the modified data. The Restore All button is now used to restore the plots to their original appearance. Then, a data point closer to the centroid of the x values, say, point 9, is selected similarly to demonstrate that such cases have comparatively smaller leverage values.
Going back to point 18 or 19, notice that when one of these is moved, the fitted line attempts to "track" its movement more closely than, say, when point 9 is moved. This can be observed either directly, by following the movement of the line while the point is being dragged, or indirectly, by noting the changes in the plot of residuals or in the parameter estimates. This effectively demonstrates that leverage measures the weight each point carries in the prediction calculations. Thus, when the leverage of an observed point is large, the corresponding prediction for that point attempts to move closer to or "track" the observed point. This forces the fitted line to pass closer to points with higher leverage.
Point 18 is deleted by first selecting it and then pushing the Delete Selection button. This will remove point 18, and straight line model is refitted to the resulting data. The change in the fit statistics is noted and the point replaced by pushing the Restore All button. Figure 6 and Figure 7 display the results of these operations. By repeating the same procedure with point 19, it will be revealed that fitted values are more sensitive to the deletion of case 19 than of case 18. Case 19 is said to be a more influential observation than case 18. This fact is reflected in the Cook's D statistics computed for the original data and shown in Figure 6; the largest value corresponds to case 19. Cook's D statistic is a measure of the influence an individual case has on the regression fit.

Examining Relationships among Regression Variables
Plotting Y against each of the explanatory variables and computing correlations among them are acceptable as simple tools for studying the relationships among these variables. For example, these may be useful in detecting collinearities involving a pair of variables. However, if incorrectly interpreted, these correlations alone may lead to misleading conclusions regarding the contribution of an explanatory variable to a regression model, particularly in the presence of other variables in the model. Moreover, in multiple regression, collinearities can involve three or more variables. To detect these relationships, simultaneous use of several diagnostic plots may be necessary.
Some of the important concepts related to understanding and interpreting relationships among regression variables that can be illustrated using the regteach4 module are: A plot of Y against an X variable, called the partial response plot, cannot be used alone to explain the contribution of the X variable to the multiple regression model.
A useful way to display the strength of the relationship between Y and an X variable in a model is to plot these against each other after removing the linear effects of the other explanatory variables from each. This plot is called the addedvariable plot or the partial regression plot.
Added-variable plots are useful in directly determining how individual cases affect the estimation of the corresponding regression parameters.
Finding the best linear fit is difficult when the predictors are highly correlated. This is reflected in high variance inflation factors (VIFs) for some parameter estimates in the fitted model. Figure 8 shows a (Y,X1,X2) spin-plot with buttons labelled Pitch, Roll, and Yaw which can be used to rotate the 3dimensional point cloud around any one of three fixed axes. Using the Yaw button, the point cloud can be rotated around the Y axis until the strongest fit to a straight line is observed on the plane of the computer screen. This corresponds to the fitting of a least squares plane by eye, and the line is the 2-dimensional view of that plane. In regteach4 a 2-or a 3-variable regression model with predetermined values for the X variables may be selected by pressing the corresponding button to illustrate concepts described above. The user has some control over how the data for the X variables are generated; the correlation between the X1 and X2 variables can be specified using the Corr slider-bar. For example, in Figure 8 the correlation has been set to 0.40 as can be observed in the scatterplot of X1 vs. X2. The symbols V, H, and O identify the variables plotted on each axis. When the 3-Var button is pushed, the user is given the choice of selecting the pair of independent variables to be plotted on the H and O axes by pushing one of the buttons X1&X2, X2&X3, or X1&X3.
The spin-plot gives the user the ability to observe various 2-dimensional projections of the data, as demonstrated in Figure 8. In particular, the more interesting projections are those that will coincide with the partial response plots. For instance, the projection shown on the spin-plot in Figure 9 is identical to the plot of Y vs. X1 shown on the the upper left corner of the scatterplot matrix. The scatterplot matrix is used in the regteach4 module to display the relevant partial response plots simultaneously. Using the Yaw button the spin-plot can be rotated to obtain another projection that is identical to the plot of Y vs. X2, i.e., the middle plot in the top row of the scatterplot matrix.   Figures 10 and 11 show examples where the bivariate plots described above have been constructed for samples in which predictors X1 and X2 were generated with correlations of 0.04 (inducing low multicollinearity) and 0.96 (inducing high multicollinearity), respectively. The partial response plots in Figure 10 indicate a strong linear relationship between Y and X2 and a weak linear relationship between Y and X1. However, the added-variable plots show that each of these relationships are strongly linear when the linear effect of the other variable is removed from the regression. Conversely, as evident from Figure  11, lack of a linear relationship in the added-variable plot does not imply that the corresponding (Y,X) variables are independent; rather it could be that in the presence of the other variables, the X variable plotted does not provide any additional explanatory power to the fitted regression model.  Figure 12 displays an example of a 3-variable model fitted to data where the variables X1 and X2 are generated to be highly correlated with each other. By examining the added-variable plots it becomes evident that X1 does not contribute additional predictive information to the regression in the presence of X2 and X3 in the model. No linear trend is apparent in the first added-variable plot in Figure 12 showing very clearly that fitting a model with both X1 and X2 in the model will cause either or both estimated coefficients to possess a large sampling variance. The VIFs for the coefficients b 1 and b 2 are relatively large indicating the usefulness of the VIF as a direct measure of the effect of multicollinearity in model estimation. The VIFs are displayed when the Regression Stat button is selected.

Transforming Nonnormality and Nonconstant Variance
Two conditions in regression are that the mean response is a linear function of unknown regression coefficients and that the errors are additive and are a random sample from a normal distribution with constant standard deviation. Graphical tools that help to check these conditions using the fitted model are generally known as residual plots. To check the normality condition one would use a normal probability plot of either the raw residuals or studentized residuals. For checking homogeneity of error variances, plots of residuals against the predicted values and residuals against each explanatory variable (X i ) are useful. When patterns in these plots indicate deviations from model conditions, one remedy often advocated is to attempt transforming (or re-expressing) the variables to better satisfy these conditions. Transforming the response variable using the Box-Cox power family of transformations is often recommended to restore normality. More traditional transformations that are designed to achieve constant error variances also promote normality in some instances.
The regteach5 module has been constructed to illustrate some of the important concepts related to the use of residual plots and transformations in regression analysis: A normal probability plot of the residuals can be used to identify features of the shape of their distribution such as skewness and whether it is long-tailed or short-tailed.
Certain patterns in the plot of residuals against the predicted values can be used to identify the form of dependence of the error variance on the mean of Y. This plot will be called the residual plot below.
Variance stabilizing transformations achieve more nearly constant error variances. These may also help restore normality to the data in some instances.
The Box-Cox power transformation can also be used to transform Y so that the transformed data may be adequately described by a normal distribution. Figure 13 shows the start-up window of the regteach5 module. In the example displayed, the Select Data Set button has been used to select a data set and the Normal Plot button employed to obtain the corresponding normal probability plot. The start-up window also contains a Simulation button that can be used to generate data simulation windows: a window to generate simulated data and a window to display the corresponding normal probability plot (see Figure 14). While in this setup, the module is said to be in the simulation mode. In this mode, random samples may be drawn from distributions with various shapes. By observing normal probability plots of these data sets, students learn to associate types of deviations from a straight-line pattern with the shapes of the underlying distributions.
Selecting an item from the Distributions pulldown menu initiates the creation of a normal probability plot of a random sample of 30 data points drawn from one of three distributions: N(0,10), Chi-squared (3 d.f.) and Student's t (3 d.f.). A plot of the density of the selected distribution overlaid with a dot plot of the actual sample data drawn is displayed in the first simulation window. New samples can be drawn by pushing the New Sample button repeatedly or by holding it down using the mouse. Figure 14 shows the examples of plots for samples drawn from each of the above 3 distributions. In the transformation mode, the regteach5 module allows the user to perform data transformations on a selected data set. After choosing a data set by pressing Select Data Set button, the user can select a variance stabilizing scheme from a pulldown menu by pressing the Var Stbl Trans button. By observing the changes made dynamically in the residual plot and the normal probability plot in response to each transformation attempted, the user will be able to select a transformation that adequately stabilizes the variance and brings the data closer to normality. Also the Power and Shift slider-bars allow the user to try Box-Cox transformations of the original data by choosing various power and shift values. The residual plot and the normal probability plot can be saved for comparison among satisfactory transformations by pushing on the Save Plot button. Figure 15 and Figure 16 show an example of each of the above plots saved for such a comparison.
In Figure 15, the residual plot indicates a dependence of the residual variance on the magnitude of the response, and the probability plot indicates some deviation from normality. The square root transformation appears to restore normality and stabilize the variance to some extent, as evinced from the bottom set of frames. In Figure 16, the residual plot shows both curvature and nonconstant variance, although the residuals do not exhibit significant deviation from normality. As shown in the bottom set of frames, a power transformation of = 0.3 appears to stabilize the variance while retaining normality.

Availability
The instructional modules and Lisp-Stat source code for the software components of our modules are available via anonymous ftp from Iowa State University. To obtain these, use the command ftp isua.iastate.edu with "anonymous.stat" as the username and "yourusername@your.email.host" as the password. This should get you into the statistics directory named "anonymous". The subdirectories "Teach", "RegTeach", and "DsnTeach" contain Readme files describing three sets of software. If you are accessing files from a Unix host, ftp the compressed archive "regteach.tar.gz" to obtain a version that can be installed on the Unix platform. Otherwise obtain the file appropriate for your platform (e.g., regteach.exe for PC/Windows and regteach.sea.hqx for older Macs). These files are binary archives that unbundle when executed. Other Readme/Install files will be found in each package after unbundling.
Pdf formatted files of the current version of this paper and instructional module lessons designed for use with the software modules are available in the subdirectory /Docs. The subdirectory /Lessons and /Figures contain the original latex files and corresponding figures used in the lesson documents.
Lisp-Stat is freely available from umnstat.stat.umn.edu or from statlib. It is recommended that Mac OS X users adapt the Unix versions of the software; instructions to do this are available in a help file in the shell archive. Our software currently runs under Version 2.1 Release 3.52 of Lisp-Stat. We plan to make this software available from other servers on the internet such as statlib and the UCLA statistics archive.

Concluding Remarks
In this article we have extended the collection of instructional modules described in Marasinghe et al. (1996). As with the modules described in Iversen and Marasinghe (2001), the ones described here are specific to regression analysis and are much more advanced, elaborate, and complex than the ones describing elementary statistical concepts. One tenet we have attempted to follow in developing these modules is to obtain ideas and feedback from those instructors involved in teaching these topics. We hope that the research presented here will lead to exploration, refinement, and dissemination of other such modules for teaching statistics interactively. We welcome comments from potential users of these modules.
We are grateful to the statistics instructors who used the earlier sets of modules in their teaching and sent us comments and words of encouragement. The current set of modules was developed for use in courses where the primary audience is undergraduate and graduate students from disciplines other than statistics. We hope that they will enable these students to obtain an improved understanding of the underlying statistical concepts.


WebotsTM: Professional Mobile Robot Simulation

Cyberbotics Ltd. develops WebotsTM, a mobile robotics simulation software that provides you with a rapid prototyping environment for modelling, programming and simulating mobile robots. The provided robot libraries enable you to transfer your control programs to several commercially available real mobile robots. WebotsTM lets you define and modify a complete mobile robotics setup, even several different robots sharing the same environment. For each object, you can define a number of properties, such as shape, color, texture, mass, friction, etc. You can equip each robot with a large number of available sensors and actuators. You can program these robots using your favorite development environment, simulate them and optionally transfer the resulting programs onto your real robots. WebotsTM has been developed in collaboration with the Swiss Federal Institute of Technology in Lausanne, thoroughly tested, well documented and continuously maintained for over 7 years. It is now the main commercial product available from Cyberbotics Ltd.

Introduction
Cyberbotics Ltd. was founded in 1998 as a spin-off company from the Swiss Federal Institute of Technology in Lausanne (EPFL). It currently employs two people and develops Webots TM : a commercial software used for mobile robotics prototyping simulation and transfer to real robots (see Fig. 1). In 1998 and 1999 Cyberbotics developed an Aibo ® simulator for Sony Ltd. Cyberbotics has now collaborations with the Biologically Inspired Robotics Group (BIRG) and the Swarm Intelligent System Research Group (SWIS) of the EPFL through the Swiss CTI technology transfer program. Webots TM runs on Windows, Linux and Mac OS X and is intended for researchers and teachers interested in mobile robotics. It is commercially available from Cyberbotics Ltd. (http://www.cyberbotics.com).

Need for Simulation
Although the final aim is real robotics, it is often very useful to perform simulations prior to investigations with real robots. This is because simulations are easier to setup, less expensive, faster and more convenient to use. Building up new robot models and setting up experiments only takes a few hours. A simulated robotics setup is less expensive than real robots and real world setups,thus allowing a better design exploration. Simulation often runs faster than real robots while all the parameters are easildisplayed on screen. Simulations make it possible to use computer expensive algorithms that would need ages to run on real robot microcontrollers, like genetic algorithms. Finally, the simulation results are transferable onto the real robots.

Features
Webots TM has a number of essential features intended to make this simulation tool both easy to use and powerful: • Models and simulates any mobile robot, including wheeled, legged and flying robots. • Includes a complete library of sensors and actuators. • Lets you program the robots in C, C++ and Java, or from third party software through TCP/IP.

Robot and world editor
A library of sensors is provided so that you can plug a sensor in your robot model and tune it individually (range, noise, response, field of view, etc.). This sensor library includes distance sensors (infra-red and ultrasonic), range finders, light sensors, touch sensors, global positioning sensor (GPS), inclinometers, compass, cameras (1D, 2D, color, black and white), receivers (radio and infra-red), position sensors for servos, incremental encoders for wheels. Similarly, an actuator library is provided. It includes differential wheel motor unit, independent wheel motors, servos (for legs, arms, etc.), LEDs, emitters (radio and infra-red) and grippers. With Webots TM , you can create complex environments for your mobile robot simulations, using advanced hardware accelerated OpenGL technologies, including lighting, smooth shading, texture mapping, fog, etc. Moreover, Webots TM allows you to import 3D models in its scene tree (see Fig. 2) from most 3D modelling software through the VRML97 standard. You can create worlds as large as you need and Webots TM will optimize them to enable fast simulations. Complex robots can be built by assembling chains of servo nodes. This allows you to easily create legged robots with several joints per leg (as shown in Fig. 3), robot arms, pan / tilt camera systems, etc. For example, you can place several cameras on the same robot to perform binocular stereo vision, or 360 degree vision systems.

Realistic Simulation
The simulation system used in Webots TM uses virtual time, making it possible to run simulations much faster than it would take on a real robot. Depending on the complexity of the setup and the power of your computer, simulations can run up to 300 times faster than the real robot when using the fast simulation mode. The basic simulation time step can be adjusted to suit your needs (precision versus speed). A step-by-step mode is available to study in detail how your robots behave Simulating complex robotic devices including articulated mechanical parts requires precise physics simulation. Webots TM relies on ODE (Open Dynamics Engine) to perform accurate physics simulation wherever it is necessar y For each component of a robot, you can specify a mass distribution matrix (or use primitives for simple geometries), static and kinematic friction coefficients, Fig. 2. Scene tree editor window  Servo devices can be controlled by your program in torque, position or velocity. The control parameters for the servo can be individually adjusted from your controller program. The graphical user interface of Webots TM allows you to easily interact with the simulation while it is running. By dragging the mouse, you can change the viewpoint position, orientation and zoom using the mouse wheel. Pressing the shift key while dragging the mouse allows you to move or rotate objects. This feature facilitates interactive testing.

Programming Interface
Programming your robot using the C language is as simple as shown in Fig. 4. In this example, the robot is a differential wheeled robot equipped with an infra-red distance sensor named "ir" looking forward. The robot will stop moving if the distance sensor detects an obstacle and restart moving when the obstacle is no longer detected. A similar Java programming interface is also included. Moreover, any Webots TM controller can be connected to a third party software program, such as MatLab®, LabView®, Lisp®, etc. through a TCP/IP interface. Research experiments often need to interact automatically with the simulation. The supervisor capability allows you to write a program responsible for supervising the experiment. Such a program can dynamically move objects, send messages to robots, record robot trajectories, add new objects or robots, etc. The supervisor capability can be used in computationally expensive simulations where a large number of robot configurations and control parameters have to be evaluated, as in genetic evolution, neural networking, machine learning, etc.

Transfer to real robots
Once tested in simulation your robot controllers can be transferred to real robots: • Khepera® and Koala®: cross-compilation of C Webots TM controllers and remote control with any programming language (see Fig. 6). • Hemisson®: Finite state automata graphical programming with remote control and autonomous execution modes. • LEGO® Mindstorms®: cross-compilation for RCX of Java Webots TM controllers based on LeJOS. • Aibo®: cross-compilation of C/C++ Webots TM controller programs based on the Open-R SDK (see Fig. 5). • Your own robot: The Webots TM user guide explains how to build your own Webots TM cross-compilation system for your very own robot.

Documentation and support
Webots TM comes with a complete documentation, including two printed manuals with color covers. This documentation is also included on the Webots TM CD-ROM in both PDF and HTML format. The Webots TM user guide explains how to install and get started with Webots TM . This manual includes a step-bystep tutorial for modelling and programming your own robot, and describes a number of sample experiments included on the CD-ROM. It explains the basic principles of Webots TM and shows you how to transfer your programs to real robots. The Webots TM reference manual contains everything you need to develop your Webots TM application. It provides a complete description of all the objects you can simulate with Webots TM , including robot bases, sensors, actuators, simple objects, etc. The programming interface is completely documented with examples. Functions are sorted by categories. Finally a number of example worlds and controllers are provided on the CD-ROM which can serve as a starting point for developing your application. Webots TM users take advantage of the Webots TM users community through a support mailing list. Most questions are answered within 24 hours by Cyberbotics support services. All Webots TM licenses include one year of personalized user support and free upgrades via the Internet.

Examples of Applications
Webots TM has been used by more than one hundred universities and research centers worldwide since 1998 for both education and research purposes. Education applications allow the students to get started with robotics, 3D modelling, programming, artificial intelligence, computer vision, artificial life, etc. with an integrated tool. It is easy to implement a virtual robot contest, like a soccer contest or a humanoid locomotion contest, based on Webots TM , which is highly motivating for the students. All the articles listed in the last section of this paper refer to different research applications of Webots TM . It would be out the scope of this paper to present each application. However, we may attempt to classify them into several major categories: • The multi-agent simulations category includes research experiments where several robots cooperate to reach a global goal (see Fig. 7). • The artificial intelligence category attempts to validate psychology hypothesis, mainly learning, by simulating intelligent mobile robot behaviors. • The control research category involves developing efficient control algorithms to perform complex mobile robot motion. • The robot design category aims at shaping mobile robots defining the position and properties of its sensors and actuators. This is especially useful to investigate new wheeled, legged or flying robots.