Clustering and PCA for Reconstructing Two Perpendicular Planes Using Ultrasonic Sensors

In this paper, the authors make use of sonar transducers to detect the corner of two orthogonal panels and they propose a strategy for accurately reconstructing the surfaces. In order to point a linear array of four sensors at the desired position, the motion of a digital motor is appropriately controlled. When the sensors are directed towards the intersection between the planes, longer times of flight are observed because of multiple reflections. All the concerned distances have to be excluded and that is why an indicator based on the output signal energy is introduced. A clustering technique allows for the partitioning of the dataset in three clusters and the indicator selects the subset containing misrepresented information. The remaining distances are corrected so as to take into consideration the directivity and they permit the plotting of two sets of points in a three-dimensional space. In order to leave out the outliers, each set is filtered by means of a confidence ellipsoid which is defined by the Principal Component Analysis (PCA). The best-fit planes are obtained based on the principal directions and the variances. Experimental tests and results are shown demonstrating the effectiveness of this new approach.


Introduction
The aim of perceiving an environment is ambitious in mechatronics because it lets a robot operate avoiding collisions. Expensive devices, such as laser sensors or camera systems, are often preferred to ultrasonic sensors in spatial reconstructions owing to the disadvantages of sound propagation. Although the latter light-based sensors are very accurate, they are not able to work in some conditions such as smoky rooms, or in the presence of light-absorbing or shining obstacles [1]. Therefore, ultrasonic sensors, also called in-air sonar sensors, are mostly used in these particular scenarios, even if their low-cost often leads researchers more and more to test their use in common home-like environments. These sound-based sensors spread mechanical waves through the air and they wait for the echoes. The elapsed time between the transmission start and the reception is the Time Of Flight (TOF) and it is proportional to the travelled distance. Several TOF measurements are generally required to estimate the distance and different methods may be compared in a statistical way [2]. The TOF values depend on the speed of sound and, as such, digital signal processing techniques prove to be a necessity for compensating for the speed variations due to temperature or other atmospheric conditions [3]. In industry, sonars are widely employed for controlling the liquid level in containers taking advantage of the sound reflection on the surface of liquid. The authors recently designed a mechatronics device for measuring the water level in a tank and for estimating the volumetric flow [4]. A lot of research is aimed at interpreting the collected data for localizing and classifying some target obstacles [5]. Other works allow for the detection of the size and shape of the reflecting object using neural networks and considering amplitude, frequency and time structure [6]. Moreover, the echoes may be processed and from this energy, duration and range can be extracted so as to characterize the roughness and the orientation of the reflecting surface [7]. In some cases, the reflector geometry is such as to make the use of a multisensory sonar system more preferable than a single sensor. These systems may be advisable for tracking an object and identifying it [8] especially in conditions in which it is not possible to use other sensors. Since the sound wave front is spherical, more sensors have to be fired to locate a spherical target and to estimate its curvature [9]. Some further studies are concerned with the conversion of geometrical shapes into quantitative values. The goal of detecting openings in a wall is achieved by means of the PCA which makes it possible to attribute a scalar value to a set of points [10]. The authors have also tackled a similar problem and they have proposed a model for evaluating the regularity of a wall and for estimating its orientation [11]. Sonar sensors may be utilized to build a map of the surrounding environment. A recent processing approach for ultrasonic arc maps is proposed in [12] and compared to six other existing traditional techniques [12]. The difficulty of equipping robots with sonars is frequently due to the wide sound beam propagated, which is responsible for TOF measurements and is significantly affected by many factors, such as specular reflections and scattering. In previous studies, the authors used a new model of ultrasonic sensors able to produce an analogue sonic signal (full waveform) as well as a logic signal including information on the TOF. The authors implemented an array of four sensors moved by a digital modular actuator [13] and then utilized it to scan and reconstruct an L-shaped obstacle [14]. After this, they carefully examined the physical phenomenon of multiple reflections in the zone close to the intersection and they introduced a new reconstruction strategy based on the Fuzzy C-Means (FCM) algorithm and on the RANdom SAmple Consensus (RANSAC) fitting [15]. They have brought together the previous analyses and observations in such a way as to define a new procedure in which each best-fit plane is obtained by statistical considerations. Numerous tests have confirmed that the energy of the analogue sonic signal increases if the sensor is pointed at the intersection and a new indicator of the relating motor position is herein introduced. This indicator proved to be preferable in more cases particularly when some sensors are not able to catch reflection. The values of this indicator, the distances and the motor positions permitted us to define a cloud of points which was partitioned into three clusters by the FCM algorithm. One of these clusters contains spurious distances due to multiple reflections and it is selected by the maximum value of the indicator. The remaining sets have been used to plot some points in a threedimensional space taking into account the directivity diagrams. The aim of excluding eventual outliers has been pursued by calculating the covariance matrix of each set. The eigenvectors of this matrix give the directions of the axes of an ellipsoid, whose centroid corresponds to the mean point and the eigenvalues allow for the fixing of the length of the axes so as to contain all the inliers in it. The new covariance matrix of each filtered set permits us to obtain the normal vector of each plane and the equation may be fully determined assuming that the plane passes by the new mean point. The aim of accurately reconstructing two orthogonal planes by using ultrasonic sensors is very interesting in industrial robotics applications. For instance, fabrication processes often consist in joining together two parts by using welding robots [16][17]. The automatic welding solutions are more and more desired to reduce imperfections and to improve the final quality [18]. Since the welding process is complex, the interaction humanrobot is generally essential. The proposed reconstruction strategy could be applied to the specific problem of joining two flat surfaces together. In this way, a welding robot could become aware of the obstacle's location and the need for human intervention could be limited. Many experiments were carried out and the most meaningful ones are herein discussed.

The Sonar Sensor Array for Estimating the Distance
The rotating array is moved by means of a servo modular actuator Dynamixel AX-12+ [19]. This smart system consists of a precision motor and a control circuitry which is connected to a computer in order to drive the rotation of a frame by using the Universal Serial Bus (USB) [19]. The actuator was placed on an extendible professional tripod and its shaft was constrained to a plastic thick bar on to which four sensors were screwed. The linear array is composed of two types of ultrasonic sensors: Hagisonic AniBat HG-M40DAI and HG-M40DNI [20]. These sonar sensors are both transmitters and receivers of ultrasonic waves at a frequency of 40 kHz [20]. The HG-M40DAI sensor spreads an anisotropic beam and it may detect obstacles up to 4 m away with a horizontal directivity of about 150° and a vertical directivity of about 60°-70° [20]. The HG-M40DNI sensor spreads an isotropic beam up to 9 m and it has a directivity of 10° [20]. The capsules are placed alternately at a distance of ds = 8 cm. The rotating device is illustrated in Figure 1, where the HG-M40DAI models are labelled with the numbers "1" and "3" while the HG-M40DNI models are labelled "2" and "4". The sensors are able to rotate around the x axis of the reference system indicated in Figure 1. Whenever an input pulse is applied to an opportune wire (several shapes are available), the ultrasonic transmission is driven at its rising edge [20]. In the proposed tests the chosen input signal is rectangular with a frequency of 12.5 Hz, a duty cycle of 4% and an amplitude of 4 V. The use of a rectangular signal is recommended by the sensor producer [20] and it could be easily obtained by means of an oscillator circuit. Moreover, preliminary calibration tests showed the sensor operates well with the chosen input signal. All the sensors are powered by a recommended voltage of 12 V and a red intermittent Light Emitting Diode (LED) indicates whether it is working correctly. When the reflected wave is detected, the sensor responds to the input pulse with a subsequent output pulse. The TOF may be estimated by the period elapsed from the rising edge of the input pulse to the rising edge of the output pulse. The sensors were simultaneously triggered so that the waves were added together and all the output signals were acquired by the acquisition board NI-USB 6210 [21]. They were suitably sampled and processed to extract the desired information. In case of direct reflection, indicated by tiθ, the calculated time by the i th sensor when the motor position is θ (angle measured from the y axis to the frame direction in the yz plane), the distance may be obtained by (1) where v is the estimated speed of sound (v ≈ 345 m/s) and h is a correction factor (h = 0.922) for considering the delay of system.
In order to obtain a reliable distance estimation, the time tiθ may be obtained as the median value of several consecutive measured values. In any case, the variation in measurements taken under the same conditions is very small and so the measurements are very repeatable [13][14][15]. This variability test may be easily verified by considering the standard deviation and observing that its value is very low. The time values are very close to the mean.

The Experimental Tests and the New Indicator Based on the Signal Energy and on the Distance
The experiments are performed by using the propagation of ultrasonic waves towards the obstacle represented in Figure 2. The target obstacle is obtained by placing some wooden panels on a metal structure.
In this way, two orthogonal flat surfaces may be reached by the large beam propagated and the regions of two orthogonal planes may be individualized. In spite of the multiple reflections, the purpose of this work is to draw the planes from the measured distances in each motor position. One of the real planes is normal to the y-axis and the other one is normal to the z-axis. The proposed tests are related to the same geometric configuration. The real plane normal to the y-axis has a minimum distance dy = 157 cm from the origin of the reference system. The other plane, which is normal to the z-axis, has a minimum distance dz = 106 cm from the same point ( Figure 2). The digital motor was programmed to make a step-by-step rotation in the angular interval from -10° to 110°. The frame position was directly obtained from the control circuitry as feedback information. Two different tests were carried out in the configuration considered here with a great number of motor steps: the first test consisted of 184 steps and the second test consisted of 220 steps. Ultrasound waves were spread toward the obstacle and they returned back to the source only after the goal position was attained and when the motor was stopped. This is essential to avoiding undesired effects concerned with the motion of the source, such as the Doppler shift. The necessary experimental time for the completion of the scanning process includes the time needed to bring the sensors to the goal position and the time needed to acquire the output signals (total acquisition period). It was about 157 s for the first test and about 187 s for the second test. These times were due to the great number of rotation steps considered and to the total acquisition period of 0.8 s. The choice of such an acquisition time allowed for a reliable distance estimation from the median of the values. For the sake of reducing the scanning time, the tests can also be carried out by requiring an acquisition period of 0.08 s (only one time value for sensor in each position). In this case the scanning time goes down to about 24 s for the first test and about 29 s for the second test. Therefore, the scanning time can be reduced based on the fields of application. The time tiθ, returned by the i th sensor for all positions in the first and second test is plotted in Figures 3 and 4 respectively. In order to take statistical advantage of all four time values obtained in position θ, if a time value differs from the median value tθ = median(t1θ, t2θ, t3θ, t4θ) by less or more than 20%, it is replaced by such a median value. As expected, many spurious time values are obtained from the zone close to the intersection. The angular position which corresponds to the intersection is θ0 = atan(dz/dy) ≈ 34° (critical position). In order to exclude all the undesired values, one could choose to fix a time threshold but such an approach does not allow us to estimate θ0 and it is not advisable in the case of an unstructured environment.
The authors intend to suggest a powerful strategy which can be generally used to reconstruct the two planes even if reflection losses or unwanted reflections occur because of irrelevant obstacles interposed between the scanner and the target. Therefore, the proposed approach could be opportunely considered in mobile robot applications. In these cases, the target may be represented by room walls on which some objects are hanging (such as pictures or lamps). In order to observe the internal crosstalk and the reflected sonic signal (analogue signal or full waveform), a wire was appropriately connected to the analogue signal test point of each sensor board [20]. The other wire extremity was plugged into the acquisition board NI-USB 6210 [21]. Let viθ(t) be the full waveform given by the i th sensor in position θ, the signal energy Eiθ (expressed in V 2 s) may be calculated by (2), where T is the total acquisition period (T = 0.8 s in the experiments proposed).
Considering the sum of the energy values for the assigned position θ, an energy level may be defined by (3), where N is the number of sensors (N = 4).
The choice of considering the sum of the energy values is justified by the fact that the energy trend is quite similar for each sensor [18]. In Figure 5 the energy level is plotted for all the positions of the second test. The energy level is high in cases of direct reflection, which happens for θ ≈ 0° and θ ≈ 90° because of the high analogue signal intensity. Moreover, the energy analysis brings out another very interesting aspect. When the beam is propagated towards the cavity due to the intersection between the panels, the signal energy increases because of constructive interference and another meaningful central peak is observed. Let us denote by dθ the median distance, which is calculated as dθ = hvtθ/2, and let us consider the function W given by (4).
This function of θ allows us to automatically select the peak concerning the critical position because it brings together both the effect of the greater energy and the effect of the greater distance measured. If θ* is the estimated value of θ0, then θ* is such that W(θ*) = 1 (greater than any other value). One should note that the normalization of the product between the energy and the squared median distance is not necessary to detect the angular position corresponding to the plane intersection. However, this adjustment facilitates the next task of partitioning the data in order to exclude the misrepresented distances. In cases of different values for the total acquisition period (for example T = 0.08 s), a similar trend for W was verified. The indicator values are plotted in Figure 6 for the first test and in Figure 7 for the second test. In both cases, the critical position is fairly well estimated (θ* ≈ 36°). This new indicator of the critical position proved to be better than the previous indicator introduced in [15] which was based on the dependence of the energy on the distance. In fact, when the scanning environment is not structured, reflections against non-target obstacles in the zone close to the intersection may entail an absence of the peak. In these situations, more sensors are not able to detect the target because the reflected wave amplitude is too low or the reflected wave is directed somewhere else (absence of the output pulse for estimating the TOF). The implemented code was developed in such a way as to assign the maximum detectable distance to the reflection loss. This means that, although no energy level peak is detected, the critical position can be nonetheless estimated because the indicator presents a peak due to the great distance considered. Some tests confirmed that reflection loss often takes place in a motor position close to the critical one. Therefore the use of the indicator W turns out to be very useful in these cases.

Reconstruction Strategy Based on Cluster Analysis and PCA
The possibility of estimating the critical angular position led the authors to use it to exclude the misrepresented distances stemming from multiple reflections. Let us consider all the points (θg, diθ, W(θ)), where θg is obtained by the indexed position which is the requested address value to define the motor instruction packet [15,19]. The aim of partitioning the cloud of points is pursued by means of the FCM clustering. This soft clustering allows the points to belong to three clusters and the cluster related to the critical position is left out. In this way, the original set is filtered and the spurious distance values are removed. The FCM algorithm was developed in [22] and generalized in [23]. This classification method comes from a variation of the Lloyd's algorithm which partitions the data space into a structure known as a Voronoi diagram [24]. The Fuzzy clustering is often used in machine learning, pattern recognition, image processing and bioinformatics. The proposed approach is based on the minimization of the Euclidean distances between the points and the cluster centroids through the iterative optimization of an objective function. The grade of membership of each point in each cluster is expressed by a matrix. This membership matrix was chosen as a random matrix at the first iteration and it was updated in each following iteration as long as the cost function of the improvement was less than a minimum threshold. This is equivalent to specifying a condition on the Frobenius norm of the difference between two consecutive membership matrices [15]. The fuzziness coefficient was fixed equal to 2 for updating and the minimum amount was imposed equal to 10 -5 . These options let us find the best location for the three sets in several experiments. The clusters of the two considered tests were mapped in the (θ, diθ, W(θ)) reference system and then they were plotted with different colours as depicted in Figures 8 and 9 respectively. The excluded cluster is the red one containing the point concerned with the critical position (illustrated as a black point).  The distances diθ of the blue and green clusters may be used to define two sets of three-dimensional points in the reference system of Figure 1. Let (xiθ j , yiθ j , ziθ j ) be the generic point of the j th set such that j = 1 if θ < θ* and j = 2 if θ > θ*, the coordinates are calculated by (5) in such a way as to consider the directivity of the wide beam propagated. In the following equations, the minimum distance from each sensor capsule to the x-axis is indicated by r. Its value is about 14.5 cm for the implemented ultrasonic scanner.
The reflection scheme is indicated in Figure 10, where the effect of multiple reflections is also represented. The j th set of the three-dimensional points allows us to arrange the data block matrix Xj which is defined by (6), where the i th block is built as in (7) by considering all the ni j relating motor positions.
Please note that all elements on the first column of Xij are equal. The covariance matrix Cj of Xj is given by (8), where mj is the number of rows of Xj calculated by (9), wj is the empirical mean vector expressed in (10) and vj is a row vector of all those with mj elements.
The covariance matrix Cj permits us to compute the PCA with a view to excluding the eventual outliers from each set and to getting useful information on the orientation of the planes. The PCA consists in an orthogonal transformation in order to analyse the data [10,11]. This method is frequently used in computer vision for image compression. Let λx j be the eigenvalue of Cj associated with the unit eigenvector ux j = (ux1 j , ux2 j , ux3 j ) of Cj whose rounded components are (1,0,0), let λy j be the eigenvalue of Cj associated with the unit eigenvector uy j = (uy1 j , uy2 j , uy3 j ) of Cj whose rounded components are (0,1,0) and, finally, let λz j be the eigenvalue of Cj associated with the unit eigenvector uz j = (uz1 j , uz2 j , uz3 j ) of Cj whose rounded components are (0,0,1). Moreover, let (x0 j , y0 j , z0 j ) be the mean point of the j th set. The exclusion of the data points which are further away from the mean helps us pursue the aim of a robust estimation. The inliers are defined as all points contained in the ellipsoid which have the parametric equations given by (11), where the matrix Rj is defined in (12), the parameter ψ varies from 0 to 2π and the parameter φ varies from 0 to π. The length of the ellipsoid axes depends on the reliability factors kx and k which were chosen in such a way as to benefit from all the sensory information (kx = 3) and to reduce the confidence along the directions of uy j and uz j (k = 1). j cos sin R sin sin The data points and the two confidence ellipsoids, which were obtained for both the tests, are represented in Figures 11 and 12. The points which are inside the closed surfaces are plotted in Figures 13 and 14.    Each excluded point of the j th set permits the removal of the corresponding row from the matrix Xj so as to obtain a new matrix Xj*. This reduced data matrix allows us to calculate a new covariance matrix Cj* in the same way as described above. The aim of drawing the planes was soon reached. In fact, the smallest eigenvalue of Cj* is associated with a characteristic vector whose direction is normal to the fitting plane. Let us suppose that the selected unit eigenvector is nj = (nx j , ny j , nz j ), the equation of the plane may be expressed by (13), where the coefficient q j is calculated as in (14) assuming that the plane passes through the median point (mx j , my j , mz j ) of the considered set.
The use of the PCA allows us to achieve the very accurate reconstructions illustrated in Figure 15 for the first test and in Figure 16 for the second test. The fit planes, depicted as coloured full surfaces, are compared with the expected ones, depicted by mesh.  To appreciate the validity of the proposed approach in comparison to the strategy recently introduced in [15], the equations of the fit planes were compared with the outcome of a RANSAC function [25] applied to the original sets of points calculated by (5). The RANSAC consists of a least squares estimator which similarly uses a restricted range of points (inliers). The resulting planes are plotted in Figure 17 and in Figure 18 for the two tests.
As well as the previous reconstructions, the outcoming planes are compared with the real locations.  The comparison between the estimated parameters is expressed in Table 1 for the first test and in Table 2 for the second test. In both experiments they are very close to the real ones. In order to obtain the inclination of the reconstructed plane with respect to the three axes, the angle between the normal vector to the estimated plane and the normal vector to the real plane was also calculated. The angle β j for each reconstruction is also indicated in the tables to demonstrate the effectiveness of the reconstruction strategy. Both RANSAC and PCA may be used to produce results which match the set of observations very well. The more motor steps were executed in an experiment, the better the final RANSAC-based reconstructions came out because of the greater number of inliers. This improvement is not so striking because of the great number of motor steps considered in the two tests. The angle β j is close to 0° in all the cases but the PCA-based approach is also able to operate better than RANSAC on smaller sets of point (such as the first test). On the other hand, the PCA for model fitting requires a little bit more computation effort due to the calculation of the parametric equations for the ellipsoids. The total execution time for the proposed approach was about 0.23 s for the first test and 0.29 s for the second test. The code was executed on a computer with 3.40 GHz CPU and 8 GB RAM. Hence, the clustering and PCA-based reconstruction strategy is preferable above all because the motor programming of many rotation steps causes a longer scanning process. The full FCM and PCA-based reconstruction strategy is explained by the following flow chart in Figure 19.

Conclusion
Nowadays, technology more and more deals with the design of robots that can help humans in several fields, such as manufacturing. A lot of research focuses on investigating new ways to let a robot become aware of its environment. Recent studies aim to endow the new generation systems with sensory abilities in such a way as to assure safety in operating. In this work the authors use an ultrasonic scanner of their own design in order to spread sound waves towards a target obstacle with the purpose of reconstructing it. When a sound beam is propagated for obstacle detection, the main difficulty stems from some undesired physical phenomena which may take place. Specular reflections may be responsible for incorrect detection. The careful examination of the related effects allows the development of a reconstruction strategy in case of multiple reflections. Since the scanning of L-shaped obstacles brings out further drawbacks to the diffuse propagation, such as undesired reflections, reflection losses and scattering, the authors take advantage of a soft clustering technique and of a statistical analysis to exclude all the spurious information and fit the planes. The proposed approach allows us to attribute a confidence to the data and each plane is drawn by the principal coordinate analysis. The objective of reconstructing the planes could be similarly attained by using laser scanners but this could be unsuitable in some scenarios. For instance, a robot could be equipped with the implemented scanner in such a way as to weld together two planes. In this case, a common laser device could fail because of the smoke. On the contrary, the sonars do not feel the effect of low visibility and they may also work in environments surrounded by high reflecting objects. The final reconstructions are very good and they provide encouragement for future studies.