Reconstruction and uniqueness of moving obstacles

09/01/2013
by   Kamen M. Lozev, et al.
0

We study the uniqueness and accuracy of the numerical solution of the problem of reconstruction of the shape and trajectory of a reflecting obstacle moving in an inhomogeneous medium from travel times, start and end points, and initial angles of ultrasonic rays reflecting at the obstacle. The speed of sound in the domain when there is no obstacle present is known and provided as an input parameter which together with the other initial data enables the algorithm to trace ray paths and find their reflection points. The reflection points determine with high-resolution the shape and trajectory of the obstacle. The method has predictable computational complexity and performance and is very efficient when it is parallelized and optimized because only a small portion of the domain is reconstructed.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 1

page 2

page 3

page 4

11/28/2011

Shape and Trajectory Tracking of Moving Obstacles

This work presents new methods and algorithms for tracking the shape and...
08/06/2019

Fast Time-optimal Avoidance of Moving Obstacles for High-Speed MAV Flight

In this work, we propose a method to efficiently compute smooth, time-op...
09/24/2018

Rectilinear Shortest Paths Among Transient Obstacles

This paper presents an optimal Θ(n n) algorithm for determining time-mi...
03/25/2022

Moving Obstacle Avoidance: a Data-Driven Risk-Aware Approach

This paper proposes a new structured method for a moving agent to predic...
06/09/2017

An Efficient Algorithm for Computing High-Quality Paths amid Polygonal Obstacles

We study a path-planning problem amid a set O of obstacles in R^2, in wh...
04/27/2021

Multifrequency inverse obstacle scattering with unknown impedance boundary conditions using recursive linearization

We consider the reconstruction of the shape and the impedance function o...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Let be a reflecting convex moving obstacle with a smooth boundary in a domain and . Consider an ultrasonic wave, or signal, described by the wave equation

(1)

where is the variable speed of sound in . We consider as positive the speed of sound at each point of the domain which is not inside an obstacle and define

and

for . When then and therefore we consider that the speed of sound in is positive and known when there is no obstacle present.

We model the signals as rays and look for solutions of the wave equation. These solutions of the wave equation [1] are called rays or ray solutions [2]. We model to be an environment without caustics.

Let in . Suppose that for all t we are given all integrals where are broken rays in reflecting at . A broken ray is a ray reflecting at the obstacle and is defined as the union of two rays and in starting at the observation boundary and intersecting at the obstacle’s boundary. Then as we know correspond to signal travel times in a medium with speed of sound . The shape and trajectory reconstruction problem is to find given the sets of travel times where , the initial and end positions and take off angles of the rays , and the speed of sound in the whole domain when there is no obstacle present.

In our model for the shape and trajectory reconstruction problem rays start at signal transmitters and end at signal receivers with known locations in the observation boundary . In addition, we model the rays to have known initial conditions: the initial zenith and azimuth angles at the transmitter as well as the times when signals are sent are recorded and are known. Receivers can record the times when signals are received and these times are known as well. The combined information from transmitters and receivers provides the data or data points

for the shape and trajectory reconstruction problem where where and are the initial incident and azimuth angles of the ray with index k from its transmitter, , , are the coordinates of the transmitter endpoint of the ray, and , , are the coordinates of the receiver endpoint of the ray, is the travel time for the signal and is a frequency of the signal.

The ray paths of unbroken rays with known initial conditions are solutions of a system of equations used in the Shooting Method for two-point seismic ray tracing[3]:

(2)
(3)
(4)
(5)
(6)

Systems of equations that present an initial value formulation for the ray equations have origins in acoustics [4] and are used in algorithms for seismic ray tracing [5, 6, 7]. In other imaging fields, such as non-destructive testing and biomedical imaging, non-linear ultrasound is studied with a focus on frequency methods[8, 9]. This work provides a mathematical definition and solution of the shape and trajectory reconstruction problem and is focused on reconstruction and uniques of moving obstacles.

In the above system of equations (x(t), y(t), z(t)) is the ray position vector,

is the incident angle of the ray direction vector with the z axis and is the azimuth angle that the projection of the ray direction vector makes with the positive x axis.

In order to reconstruct the obstacle, we consider the speed of sound to be positive and known throughout when there is no obstacle present and trace rays from transmitters and receivers as if there is no obstacle. When the sum of ray travel times at an intersection point of a transmitter and a receiver ray is equal to the travel time from the corresponding data point, we infer that the intersection point could be a reflection point from . The algorithm is described in the paper on shape and trajectory reconstruction of moving obstacles[10] and its operation is shown in Figure 1.

Figure 1: In order to find a solution point P that satisfies the initial conditions for a data point, the algorithm traces a ray from transmitter L with initial angles from the data point and traces rays from receiver S with all possible initial angles. When a receiver ray intersects the transmitter ray at a point P, the algorithm checks whether the sum of travel times for LP and SP equals the travel time from the data. When the sum of travel times is equal to the travel time from the data point then P is included in the solution set. This paper describes a new filtering procedure for the solution set which ensures that P is the unique reflection point for the measurement ray for the data point.

This work extends the above algorithm with a filtering phase which ensures that the reconstructed solution is unique and contains only points from . This work also extends the first phase of the reconstruction algorithm for computing all reflection points that meet the initial conditions with an adaptive computation of the time step for tracing transmitter and receiver branches of broken rays. We then analyse the conditions for existence and uniqueness of the solution.

2 Reconstruction algorithms

The reconstruction algorithms work in two phases. The first phase finds all points in that are intersection points of transmitter and receiver rays with initial conditions from the data and for which the sum of travel times of the traced rays from transmitter and receiver endpoints to the intersection point is equal to the travel time from the corresponding

data point. This section presents the first phase of the algorithm. The next section describes and analyses the second filtering phase.

The input to the following algorithm is the speed of sound for the domain when there is no obstacle present and a set of data points or ray coordinates corresponding to the initial conditions of broken rays and their travel times. The output is a set of points in reconstructed from the input data.

0:  Set of broken ray data points
0:  Speed of sound for domain

when there is no obstacle {Algorithm for Shape and Trajectory Reconstruction of Moving Obstacles} {Estimated time complexity is

where T is the number of discretization points for the time of flight, and A is the number of discretization points for the angle space}
  for all data points  do
      set this initial position to be position of transmitter
      set this initial position to be position of receiver
     
     
     
     
     repeat
        {Compute the next point on the ray from the transmitter by Runge-Kutta step and the ray tracing system 2} { is a constant value for initializing the time steps}
        
        repeat
           )
           )
           )
           )
           )
           )
           )
           )
           )
           )
           if  ( or or or or )  then
              
              done=false
           else
              done=true
           end if
        until done or ()
        if  then
           {step size is too small. exit computation for this data point and continue with next data point .}
        end if
        
        if  then
           {We are over the travel time budget . Continue with next data point .}
        end if
         point on solution of ray tracing equations with initial values for transmitter that is at time away from the transmitter L
        if !(then
           There must be a measurment error. Continue with next data point
        end if
        for all initial angles in discretized angle space of the receiver do
           
           repeat
              {Compute the next point on the ray from the receiver by Runge-Kutta step and the ray tracing system 2}
              repeat
                 )
                 )
                 )
                 )
                 )
                 )
                 )
                 )
                 )
                 )
                 if  ( or or or or )  then
                    
                    done=false
                 else
                    done=true
                 end if
              until done or ()
              if  then
                 {step size is too small. exit computation for this data point and continue with next data point.}
              end if
              
               point on solution of ray tracing equations with initial angles and and initial position S, that is time away from S
              if !(then
                 Exit this for loop and continue with next pair of initial angles from outer for loop
              end if
              if  and  then
                  {Solution for current data point found. Add to list of all solutions for all data points. Continue with the outer transmitter loop to look for more solutions for .}
              end if
              if  then
                 {We are over the travel time budget . Continue looking for a solution with the next set of initial angles .}
              end if
           until threshold for maximum number of time steps from receiver
        end for
     until threshold for maximum number of time steps from transmitter
     {No solution found for due to measurement or other errors. Continue with next data point.}
  end for

The algorithm uses a Runge-Kutta method for the RK step and is flexible to work with other time-dependent numerical methods. When the algorithm is parallelized and caching and other optimization techniques are used then its computational complexity is where T is the number of discretization points for the travel time of a broken ray. For input from one sampling time interval , the algorithm reconstructs the shape of the obstacle during this sampling interval and the trajectory of the obstacle is reconstructed when the algorithm is run on the data points for each of the sampling intervals. Resolution of the reconstruction can be very high because the reconstruction method allows collection and processing of a large number of data points corresponding to different points from . Reconstruction with high resolution by the above algorithm of a neighborhood of points from a moving obstacle is shown in Figures 2, 3, 4, 5 and 6.

Figure 2: Reconstruction of a line segment in an environment with speed of sound . Both branches of each broken ray are curves.
Figure 3: Reconstruction of the same line segment as in Figure 2 when moving to a new location in the same environment with speed of sound . Both branches of each broken ray are curves.

3 Existence and uniqueness of a reconstructed point

For a transmitter at and receiver at the travel time between L and S along a broken ray with segments and and reflection point is

(7)

For a fixed and constant speed of sound c the above equation implies that the set of points P is an ellipsoid with focci and . This can be seen by multiplying both sides of the equation by which leads to

(8)

which is the equation of an ellipsoid with focci and .

Therefore, for constant speed of sound and a data point with transmitter , receiver , initial angles and and travel time t, a unique point P can be reconstructed because a ray from L with initial angles and travelling along a straight line, because of the constant speed of sound, will intersect the ellipsoid in exactly one point P.

For variable speed of sound the set of points is a surface which is not necessarily convex. In this case a ray from L with given initial conditions and travelling along a curve could intersect the surface in more than one point. In other words there can be several rays from the receiver that intersect the transmitter ray at different points such that for each intersection point the sum of travel times from transmitter and receiver to the intersection point is equal to the total travel time.

In order to find the unique reflection point for the measurement ray, the above algorithm is extended by adding a second filtering phase after the first phase of the algorithm. The first phase finds all solution points for which the sum of travel times from transmitter and receiver to the solution point is equal to the travel time from the data point. The filtering phase finds the unique reflection point by reconstructing the shape via several pairs of transmitters and receivers and selecting those points which have been recontructed by a sufficiently large threshold number , , of (transmitter, receiver) pairs. This implies that the data must contain at least q data points with different transmitter receiver pairs for each reconstructed point. Therefore, the observation boundary must contain a sufficient number of transmitters and receivers that are located so that each point in is seen from at least q transmitter receiver pairs. The algorithm with the combined first and second phases is as follows.

0:  Set of broken ray data points
0:  Speed of sound for domain when there is no obstacle {Algorithm for Shape and Trajectory Reconstruction of Moving Obstacles} {Estimated time complexity is where T is the number of discretization points for the time of flight, and A is the number of discretization points for the angle space and N is the number of data points. It is possible to implement the filtering phase with computational complexity which leads to computational complexity of for the whole algorithm.}
  for all data points  do
      set this initial position to be position of transmitter
      set this initial position to be position of receiver
     
     
     
     
     repeat
        {Compute the next point on the ray from the transmitter by Runge-Kutta step and the ray tracing system 2} { is a constant value for initializing the time steps}
        
        repeat
           )
           )
           )
           )
           )
           )
           )
           )
           )
           )
           if  ( or or or or )  then
              
              done=false
           else
              done=true
           end if
        until done or ()
        if  then
           {step size is too small. exit computation for this data point and continue with next data point .}
        end if
        
        if  then
           {We are over the travel time budget . Continue with next data point .}
        end if
         point on solution of ray tracing equations with initial values for transmitter that is at time away from the transmitter L
        if !(then
           There must be a measurment error. Continue with next data point
        end if
        for all initial angles in discretized angle space of the receiver do
           
           repeat
              {Compute the next point on the ray from the receiver by Runge-Kutta step and the ray tracing system 2}
              repeat
                 )
                 )
                 )
                 )
                 )
                 )
                 )
                 )
                 )
                 )
                 if  ( or or or or )  then
                    
                    done=false
                 else
                    done=true
                 end if
              until done or ()
              if  then
                 {step size is too small. exit computation for this data point and continue with next data point.}
              end if
              
               point on solution of ray tracing equations with initial angles and and initial position S, that is time away from S
              if !(then
                 Exit this for loop and continue with next pair of initial angles from outer for loop
              end if
              if  and  then
                  {Solution for current data point found. Add to list of all solutions for all data points. Continue with the outer transmitter loop to look for more solutions for .}
              end if
              if  then
                 {We are over the travel time budget . Continue looking for a solution with the next set of initial angles .}
              end if
           until threshold for maximum number of time steps from receiver
        end for
     until threshold for maximum number of time steps from transmitter
     {No solution found for due to measurement or other errors. Continue with next data point.}
  end for
  {Filter solution set.}
  for all solution points P do
     {count how many times P is reconstructed by different transmitter and receiver pairs by checking the distance between P and all other points Q in the solution set:}
     if  and the data points for P and Q have different (transmitter,receiver) pairs then
        P.count++
        remove Q from solution set because it is already counted
     end if
  end for
  for all solutions P in the filtered solution set do
     if  then
        remove P from solution set because it is an intangible solution point and not a reflection point i.e. reflection points are reconstructed by a sufficiently large number of transmitter receiver pairs.
     end if
  end for

The proof of the correctness of the filtering phase and uniquess conditions on the input data are as follows.

Theorem 3.1.

Uniqueness of a point reconstructed from multipile measurements. Each point can be reconstructed uniquely when the set of data points for every sampling interval contains at least q measurements of from q different transmitter receiver pairs where q is a sufficiently large threshold number and .

Sketch of proof: For each data point the first phase of the above reconstruction algorithm finds one or more solution points ,…,. Only one of these solution points is the unique reflection point P for ’s measurement ray. The remaining solution points for

will be referred to as intangible solution points. The conditions of the theorem guarantee that for P there are at least q data points with different transmitter receiver pairs and this implies that P will be counted at least q times by the algorithm. The probability

that any one of the intangible solution points is also an intangible solution point for another data point and its measurement ray is less than 1 and depends on the discretization of the numerical solution. Therefore, the probability that any one of the intangible solution points is counted at least q times by the second phase of the algorithm and each of these times it is an intangible solution point for the corresponding data point is less than or equal to . For sufficiently large q this probability tends to 0, therefore, with probability one, solution points that are not unique reflection points for at least one measurement ray will be filtered out by the algorithm. Therefore, the conditions of the theorem guarantee that for each a unique reflection point P with count greater than or equal to q exists because each point from the obstacle’s boundary is measured from at least q different transmitter receiver pairs.

4 Reconstruction tests

Consider a circular reflecting obstacle with Lambertian reflectance in the plane xy moving away from the origin along the line in a medium with variable speed of sound . We place a transmitter and a receiver at the origin. In this case, the domain is a circle of sufficiently large radius that contains the origin. Table 1 shows the computation by the algorithm from section 2 of the trajectory of a point on the obstacle on the line corresponding to data with different travel times from different sampling periods.

xl yl zl xr yr zr T xp yp zp
0.00 0.00 0.00 0.00 0.00 0.00 1.57 0.79 0.25 0.09 0.09 0.00
0.00 0.00 0.00 0.00 0.00 0.00 1.57 0.79 0.5 0.21 0.21 0.00
0.00 0.00 0.00 0.00 0.00 0.00 1.57 0.79 0.75 0.35 0.35 0.00
0.00 0.00 0.00 0.00 0.00 0.00 1.57 0.79 1.0 0.51 0.51 0.00
0.00 0.00 0.00 0.00 0.00 0.00 1.57 0.79 1.25 0.71 0.71 0.00
0.00 0.00 0.00 0.00 0.00 0.00 1.57 0.79 1.5 0.94 0.94 0.00
0.00 0.00 0.00 0.00 0.00 0.00 1.57 0.79 1.75 1.22 1.23 0.00
0.00 0.00 0.00 0.00 0.00 0.00 1.57 0.79 2.00 1.55 1.55 0.00
Table 1: Reconstruction of a point moving with speed on the line for a fixed signal frequency . The initial transmission angles are and . The reconstruction data for each row is from a different sampling period.

We check whether the computation of the above table by the algorithm from 2 is correct as follows. The time for the ray to reach to obstacle can be computed by the formula

Therefore,

By symmetry, for this particular example, this time t is half of the total travel time T. Then for a travel time , or , we compute . This result matches the corresponding result for xp and yp from Table 1 obtained by numerical integration. For the numerical computation gives while by the above formula . In this case, the relative error between the values computed by the algorithm and the value computed by the formula is less than 1 percent.

In order to reconstruct more points from we can vary the transmission angles. Table 2 shows that by varying , the initial angle at which we transmit rays from transmitter, we can reconstruct points on the boundary of the circular obstacle. In contrast, in Table 1 both and are constant. In order to reconstruct the boundary with higher resolution we can change the initial angles and in smaller steps. Figure 5 and Figure 6 show how changing the initial angle at the transmitter in smaller steps leads to higher resolution of the reconstruction.

The speed of the obstacle must be sufficiently slow compared to the speed of the signals during every sampling period so that for the signal travel times T from one sampling period, i.e. the times in column T where the sampling period is the same, for time period , where

(9)

the obstacle does not move by a noticeable amount. denotes the index or unique id of a sampling time period and the duration of this time period. Therefore, the duration of the sampling time period , during which rays are sent from the transmitter in order to reconstruct the boundary of the obstacle when it is approximately stationary, must be less than or equal to or

(10)

Combining the images of the obstacle from successive sampling intervals reconstructs the shape and trajectory of the moving obstacle.

xl yl zl xr yr zr T xp yp zp
0.00 0.00 0.00 0.00 0.00 0.00 1.57 0.78 1.55 1.00 0.99 0.00
0.00 0.00 0.00 0.00 0.00 0.00 1.57 0.75 1.56 1.07 0.92 0.00
0.00 0.00 0.00 0.00 0.00 0.00 1.57 0.71 1.58 1.15 0.86 0.00
0.00 0.00 0.00 0.00 0.00 0.00 1.57 0.68 1.61 1.24 0.81 0.00
0.00 0.00 0.00 0.00 0.00 0.00 1.57 0.66 1.65 1.32 0.77 0.00
0.00 0.00 0.00 0.00 0.00 0.00 1.57 0.64 1.70 1.42 0.74 0.00
Table 2: Reconstruction of points on the boundary of an obstacle in an environment with speed of sound . The initial transmission angle and the initial transmission angle varies around . The reconstruction data for all rows is from one sampling time period during which the obsacle is approximately stationary.
Figure 4: Reconstruction of points from the circular obstacle with a transmitter and receiver located at the origin in the same environment with speed of sound . The obstacle’s reflectance is Lambertian and the detected reflected ray travels to the receiver along the same path as the incident ray.
Figure 5: Reconstruction of points from the circular obstacle with a transmitter and receiver located at the origin and another transmitter and receiver located at and in the same environment with speed of sound . The obstacle’s reflectance is Lambertian and we show the reflected rays intercepting the receiver.
Figure 6: Reconstruction of portions of the boundary of the circular obstacle with a transmitter and receiver located at the origin and another transmitter and receiver located at and in the same environment with speed of sound . In constrast to Figure 5 more rays with initial angles that are closer to each other are used for reconstruction and as a result more points are reconstructed that are closer to each other i.e. the same portions of the boundary are reconstructed with higher resolution. The obstacle’s reflectance is Lambertian and we show the reflected rays intercepting the receiver.

The reconstruction tests show that for high resolution and performance it is essential that the method is run adaptively. In addition to changing the time step, the reconstruction accuracy can be tuned by changing the number of tested angles at the receiver.

0:  Set of broken ray data points
0:  Speed of sound for domain when there is no obstacle {Algorithm for Shape and Trajectory Reconstruction of Moving Obstacles} {Estimated time complexity is where T is the number of discretization points for the time of flight, and A is the number of discretization points for the angle space}
  for all data points  do
     run in parallel the algorithm from section 2 in order to reconstruct the solution points for data
     if no solution found for  then
        Run the algorithm from section 2 with finer discretization of the angle space i.e. double the number of tested angles at the receiver in order to reconstruct the point for data . Repeat until a solution is found or the angle steps become smaller than a threshold.
     end if
  end for
  Filter the solution set by the algorithm from section 3

The above algorithm adapts the angle resolution and via the algorithm from section 2 it adapts the time step.

5 Error Analysis

In floating point arithmetic

(11)

where is the time for the transmitter segment, is the time for the receiver segment and t the total travel time for the broken ray. The check from the algorithm

  if  and  then
      {Solution for current data point found. Continue with next data point }
  end if

implies that there are many points that are sufficiently close to a solution. The error is determined by the constants and and it is necessary to choose sufficiently small constants for reconstruction with high accuracy and resolution. Therefore, reconstruction of a point is unique within a ball of sufficiently small diameter which depends on the constants and .

The error of the solution is also determined by the sum of the two errors from the numerical integrations for the rays from transmitter and receiver. In addition, the error of the solution is determined by discretization errors. The initial angles at the receiver from the numerical solution belong to a finite set of initial angles that are tried. The differences between the initial angles from the numerical solution and real angles for the ray at the receiver are

(12)

and

(13)

These differences are guaranteed to be sufficiently small when the initial angles space is tested in sufficiently small equal steps.

6 Performance Optimizations

One performance optimization of the algorithm for shape and trajectory reconstruction of moving obstacles is the use of a data structure or a database for looking up points on the receiver rays for a given discretization of the initial angles space. Consider a cover of by a finite number of cubes. Let M be a cube, , centered at the origin and with sides parallel to the xy, yz and xz planes. Let be the length of one side of and divide into cubes with side length where is a natural number that determines the resolution of the mesh. Assign each of the cubes with side b from the resulting mesh a unique natural number from 1 to . All precomputed points on receiver rays are stored in a database table RT with the following schema:

pointid rayid receiverid region t x y z
Table 3: Reconstruction database table RT schema: pointid is a unique identifier for each point on a receiver ray, region is the number of the cube from the mesh that contains the point, t is the time to reach the point from the receiver, and (x,y,z) the coordinates of the point.

The region column corresponds to the number of the cube from the mesh to which point belongs. The region or cube number of a point can be defined by the function

(14)

The optimized first phase of the algorithm with a fixed time step and lookup of cached receiver rays is as follows.

0:  Set of broken ray data points
0:  Speed of sound for domain when there is no obstacle. {Algorithm for Shape and Trajectory Reconstruction of Moving Obstacles} {Estimated time complexity is when using memory databases/data structures of precomputed receiver points. T is the number of discretization points for the time of flight}
  for all data points  do
     
      set this initial position to be position of transmitter
      set this initial position to be position of receiver
     
     
     
     
     repeat
        {Compute the next point on the ray from the transmitter by Runge-Kutta step and the ray tracing system 2}
        )
        )
        )
        )
        )
        
         point on solution of ray tracing equations with initial values for transmitter that is at time away from the transmitter L
        if !(then
           There must be a measurment error. Continue with next data point
        end if
        
        for all points in and adjacent regions do
           if