DETECTION OF SPECIFIC AREAS AND DENSITIES FOR ULTRASOUND TOMOGRAPHY

In this paper we consider the problem of ultrasound tomography. Recently, an increased interest in ultrasound tomography has been caused by non-invasiveness of the method and increased detection accuracy (as compared to radiation tomography), and also ultrasound tomography does not put at risk human health. We study possibil-ities of detection of speciﬁc areas and determining their density using ultrasound tomography data. The process of image reconstruction based on ultrasound data is computationally complex and time consuming. It con-tains the following parts: calculation of the time-of-ﬂight (TOF) of a signal, detection of speciﬁc areas, calculation of density of speciﬁc areas. The calculation of the arrival time of a signal is a very important part, because the er-rors in the calculation of quantities strongly inﬂuence the total problem solution. We offer ultrasound imaging reconstruction technology that can be easily parallelized. The whole process is described: from extracting the arrival times of signals raw data feeding from physical receivers to obtaining the desired results.


Introduction
Ultrasound imaging has attracted increased interest due to the significant progress of scanning and data processing tools. In particular, in diagnostic medical procedures, including neurosurgery, in screening focal changes [Miller et al., 2012;Matthews et al., 2017]. The standard in breast cancer detection is mammography, but the use of ultrasound tomography is steadily increasing. The reasons for this are several essential factors. First of all, ultrasound examination is completely noninvasive and safe for the patient. Secondly, the proportion of women with heterogeneous or extremely dense breasts varies from 32.5% to 45.7% between regions, and the probability of an error of the traditional methods of breast screening in such cases is high [Kolb et al., 2002]. In this article, we describe the continuation of our previous research [Erofeeva et al., 2018] that is already being conducted on real experimental data.

General View of the Process
Ultrasound tomography uses various imaging techniques: image reconstructed using the reflected signal, using the pass-through signal (sound speed imaging), using the signal attenuation, and combinations of these.
We consider the through and reflected signals. The time of arrival of signals (through and reflected) is used in image reconstruction. In dense tissue the speed of propagation of the signal is higher than in the surrounding tissue. Therefore, the signal passes through such an area faster than through a homogeneous one (without specific areas of dense tissue). On the basis of such "deviations" in the time of arrival of signals, image reconstruction is later implemented. A general view of the tomographic process can be represented by the following sequence of actions: signal emission, collecting observations from sensors, extraction of signal arrival time, image reconstruction. The first stage is performed on a physical device. The device is a ring transducer with N elements (emitting sensors) (Fig. 1). Each element is an ultrasound emitter and receiver. The image is obtained by moving this ring in the vertical direction. The result of such scanning procedure are transverse "slices" of the tissues from which a 3D image is formed. When one element emits a signal, all other elements receive signals, thus the reconstruction of one "slice" accounts for N 2 signals. That is, each of the sensors in turn emits a signal and all sensors receive its signal for a fixed time. In the third stage, indications on a specific sensor at discrete points of time are used to calculate the time of arrival of a signal (Time Of Flight, TOF) to the sensor. Then TOFs are used to reconstruct the image in the fourth stage [Li et al., 2009a]. At the fourth stage, the inverse problem of image reconstruction is solved. The purpose of the inverse problem is to estimate the velocity distribution, which corresponds to the trajectories of the shortest signal passage in the region of interest. The studied area is discretized using a grid overlay.

Calculation of the Time-Of-Flight of a Signal
The differences in the speed of passage through different environments are quite small, so it is very important to extract the TOF most accurately. Considered ways to calculate TOF are: setting thresholds, comparison of a signal with a reference wave and using the Akaike Information Criterion.

Setting Thresholds
In this approach the point of arrival of the signal is considered as the moment when the signal level was above a predetermined threshold, usually dependent on the noise level. This option is poorly suited when the data have a low signal-to-noise ratio, that is, the level of noise and signal are comparable [Li et al., 2009b]. It was decided to abandon this approach due to the specifics of the data arising from the geometry of the device (ring). Since for different emitter-receiver pairs the signal level is different. For receivers-sensors that are close to the signal source, the signal level is about a hundred times lower than for diametrically opposite receivers from the signal source, that is explained by the physics of signal propagation (Fig. 2).

Comparison of a Signal with a Reference Wave
This method assumes that the signal is quite similar to a certain standard (known in advance). Comparison of the signal with the reference signal is made. Also, a measure of signal similarity is preselected. The TOF of a signal is defined as the point at which the maximum of this measure is reached [Li et al., 2009b]. Since the signals differ for different pairs of sensors and different media, this method also does not look applicable to the problem.

Akaike Information Criterion
We use Akaike Information Criterion (AIC) in the part of acoustic signals TOFs extraction based on the assumption that a wave within a certain time window can be divided into two segments -before and after the signal. This method was chosen as the best, and further work was related to its software implementation.
A modified formula for calculating the AIC was proposed in [Zhang et al., 2003] and is used for signals of different nature, including ultrasonic signals [St-Onge, 2011].
Input of the method is a time window, within which a TOF search will be conducted. Within this window, for each k, for k = 1, . . . , T , where T is the number of observations in the window, the AIC is calculated using the following formula: S(1, k) and S(k + 1, T ) --are segments obtained by splitting the window at the point k. The var(·) function has the following form: The minimum point of AIC is selected as the desired TOF.

Process Description
Pre-processing of data obtained using an ultrasound device includes the steps of: 1. Recording data for an emitter-receiver pair. 2. Search for a time window containing the true point of signal TOF. 3. Calculation of the AIC within the selected window. 4. Selection of the minimum point of the AIC as the moment of a signal arrival.

Time Window Selection
The time window within which the AIC is calculated should include the TOF of the signal. During the work various options for window selection were considered: 1. The neighborhood of the point at which the signal level on the sensor for the first time reached 2 3 of the maximum value. 2. Based on the relative position of the emitting and receiving sensors and the speed of sound in a medium quite similar to the one under study (for example, in water, since the speed of sound in water is close enough to the speed of sound in breast tissue) [Li et al., 2009b]. 3. Using the signal dispersion, calculated for smaller intervals.
Option 1 was used at the beginning of work as the easiest. The disadvantage of this approach is that if there is a stronger than the direct signal (for example, the reflected signal), the TOF will be extracted incorrectly.
Option 2 was used further as more accurate to determine the TOF of direct signal. This approach requires pre-calculated arrival times of the signal in the medium without obstacles, which causes additional difficulties.
The first two options are not suitable if not only direct signals, but also reflected ones are of interest. Option 3 was proposed to detect TOFs of all recorded signals and is as follows: − The interval width T is selected (In practice, the value T = 60 discrete times was used). − Dispersions are calculated for all intervals [x, x+ T 2 ] within the whole interval in which signal level was recorded. − The peaks of the obtained dispersion values are found, they are used as time window centers for calculating AIC and searching for TOF.

Results of Extraction TOF
The third method calculates the TOF of the signal. The Fig. 3 shows the result of this approach: the TOFs of all signals are well detected.  Table 1 shows the running time of the software implementation of TOF search. The total TOF calculation time for an emitter-receiver pair is about 60 ms. Thus, for a device with N = 100 sensors, the whole process will take about 10 minutes. This approach allows to extract all recorded signals (both direct and reflected). And also this process is very well parallelized and can be effectively implemented on GPU.

Detection of Specific Areas
Most of the existing works on similar topics are based on solving systems of linear equations [Sandhu et al., 2016;Krueger et al., 1996]. These approaches require both large computational resources and high memory costs. We are exploring a different approach that can be easily parallelized.

Description of the Image Recovery Algorithm
The algorithm is based on the fact that the time of passage of a signal between two sensors that has passed through a particular area is different from the time of passage of a signal between the same sensors in an experiment in which there are no objects to be restored. The algorithm is a sequence of the following steps: 1. The emitting sensor is fixed from all not yet considered sensors. 2. The segment between the emitter-receiver pair is painted over if the TOF of this pair is equal to TOF in the experiment without objects (that is, there were no specific areas on the signal path) (Fig. 4). 3. Steps 1 and 2 are performed for all emitting sensors not yet considered (Fig.5).

Results of Specific Areas Detection
We used several ways to implement the fill. Firstly, for the efficiency of the algorithm, a specific area was painted over. To increment the pixel values (if the reference TOF do not coincide with the TOF of this experiment), due to its speed, the Brezenham integer algorithm was used. The value of pixels that are closer to the emitting sensor is significantly higher than the value of pixels that are closer to the receiving sensor, because of the intersection of the segments. As a result, the values in some pixels were incremented more times than in others (the values are not evenly distributed). Two ways were proposed to solve this problem: with the help of additional memory for each emitter and with the help of a linear dependence of the increase in the pixel value of the distance of the emitting sensor and the receiver. At the end of the work of these implementations of the algorithm, only those pixels in which the value exceeds the threshold value are colored. To test the approach, the signal TOF were modeled. The result of the algorithm of specific areas detection is shown in the Fig. 6. The figure shows the results of the algorithm at different filling thresholds, and also used data from 50 emitters. This suggests that we can use a relatively small fraction of all the data (If sensors, for example, a couple of thousand).

Calculation of Density of Specific Areas
After the detection of specific areas, the shape and location of the object became known. The next task is to find the density of this object. The density of particular regions is related to the speed of sound in them. In the future we will look for exactly the speed of propagation in them.
There are N sensors (N = 2048) at the boundary of a certain region and M objects inside the region. Through the Snell's law on the refraction of a wave at the boundary of two medium, we want to find the speeds of sound in particular areas indicated by the blue shapes in the Fig.7. We assume that the speed of sound propagation in the medium c 0 . And in each of the regions the speed of sound propagation in the medium is c 1 , . . . , c M , moreover, all these speeds lie in a certain small, previously known range. We consider that we have a narrowly focused signal and we are able to register the number of the receiving sensor and the signal propagation time (TOF) for almost all emitting sensors. Thus, we know the pairs of interrelated points A, B and the propagation time t ab from A to B (Fig.7).
The input of the algorithm is an abbreviated table T T (dimension is N × 3) instead of the general table TOF, where dimension is N × N (reduction of 1000 times). The columns of the table T T are n i -the number of the emitter; n j -receiver number; t ij -is the TOF from i to j. It is necessary to fill in the speed map SM 2000x2000 (corresponds to the grid). In this map SM , zones are highlighted (borders are known in advance). These zones are the interior of the general field of sensors and M objects in it.
Thus, the optimization problem is posed as follows: Find such values of speeds c 1 , . . . , c M so that the calculated values oft ij (c 1 , . . . , c M ) are minimally different from the tabulated set values t ij .
We introduce the functional The calculation of values is carried out using Snell's formulas. The solution is based on an exhaustive search method for small value M and consists of the following steps: 1. Fill the sensors numbers (emitter, receiver) and TOF in the table T T . 2. In the matrix SM , set the speed 0 outside the considered region limited by sensors. 3. Set the initial approximationĉ j = c 0 for each j.
4. Fill the contours of objects in speed map SM and fill the area inside the sensors with the values c 0 . 5. The next steps are performed in a loop. Redefine the values ofĉ j for the interior of the j-th object with the values obtained from the previous step. 6. Calculate F . If the value of F is less than F min (F min determined by calculated earlier), then set the current combinationĉ j as the current estimate of the optimal solution. 7. Set the values ofĉ j so that iterates over all possible values. 8. Go to step 5, if not all values are considered.

Finding the Point of Signal TOF
To implement this algorithm for finding the velocity, you need to find interrelated points A, B.
Assumption 1. Between points A, B there is a relationship.
Assumption 2. The point of arrival of the through signal B has a maximum amplitude among neighboring receivers.
Consider an experiment in which there are no specific areas, that is, the experiment was conducted in a homogeneous medium (Fig.8). The right picture shows a graph of the dependence of signal TOF on the number of the sensor receiving the signal. And on the left graph shows the values of the maximum amplitude of these received signals. From the graphs it can be seen that the amplitudes are distributed relative to the maximum amplitude, that is, they are distributed relative to the diametrically opposite point from the point of the signal emission.
Consider an experiment in an inhomogeneous medium presented in the Fig. 9. In the study area there is an object indicated by a blue circle. Sensor number 513 emits a signal.
Then we will conduct an experiment with a specific area. The graphs of the experiment are shown in the Fig.  10. Graphics have changed: we received a drawdown on the arrival time of the signal after the signal crossed the two medium. Also the same corresponds to the graph of   the maximum amplitudes. But you can see that our local maximum is preserved.
Thus, we offer an algorithm for finding points A, B.
Compare the reference TOF graph of the experiment in water with the TOF graph of the experiment with the object and find the corridor, as shown in the Fig. 11. Then find the local maximum in this corridor, as shown in the Fig. 11 on the amplitudes graph.

Detecting the Angle of Reflection and the Point of Ray Intersection with Object
To iterate over the densities, we need to know the angle α of reflection, and the second point M of signal reflection (Fig. 7). We can find it using image-processing methods.
In the first step, the Canny edge detector is applied to the binary image. Then, detected edges of image are smoothed by Gaussian blur to obtain a solid border of the object. In the next step, the direction of the normal n is calculated. To find it, a convolution of the image with Sobel filter was used. Finally, the required angle α is calculated between the ray and the normal which are given by the points.
To detect the point M , pixel intensity values is obtained along a line that given by first point of reflection and angle of the refraction. Using this values, point M is taken as an average point between the first and last intersection of the line with the object border.

Results of Specific Areas Density
Calculation For the real experiment with N = 2048 sensors (Fig.  9), densities of the specific areas were found and are presented in the table 2. From the tables it can be seen that the obtained densities of particular regions are approximately equal for different emitters and receivers.

Conclusion
This article proposes to detect specific areas and find their density with ultrasound tomography. The whole process is described: from extracting the arrival time of a signal from raw data from a physical device to obtaining the desired results. For a series of experiments with real acoustic tomography data, calculations were performed using the proposed algorithms.