Method 1 – Traversing the Grid

For every point, we can calculate the total distance required to travel from the source to that point and then to the receiver using the distance formula on the (X,Y) coordinates of the point, the source, and the receiver.

We may now divide by velocity to get the time in seconds, and then we may again divide by the sample period to get an index for the time series vector that corresponds to this point.

Note that this value for time will not be an integer, so we must interpolate using the time series indices above and below it.

A simple linear interpolation method will give us an appropriate value that is weighted based on how close t is to t-minus and t-plus.

We finish by applying these equations to every point on the grid. This will trace out the ellipse patterns we desire for a given source-receiver pair.

Method 2 – Traverse the time vector

Taking one source-receiver pair's time vector, we first find the distance between the source and receiver. This distance is the minimum distance a signal must traverse. Since each sample in a time vector is 4ms, and we know that a wave travels at 1500m/s, the nth sample in a time vector takes n*4ms to travel n*0.004s*1500m/s. If the nth sample distance less than the distance between the source and receiver, it is ignored since it is bad data. Usually these samples have a received signal value of 0 anyway.

If the sample qualifies as valid, we must find the ellipse which satisfies the condition that the sample distance equals the distance to the reflection surface and back. Using ellipse properties, we can see that the points of this ellipse can be found relative either the source or receiver.The equation for the ellipse in polar coordinates (r,phi) from one focus is:

where e is the eccentricity of the ellipse, a is the semi-major axes, and c is the distance of the focus to the center of the ellipse. The geometry for finding a and c is shown below:

Figure 1: Ellipse geometry.

We know c as the distance between the source and receiver divided by 2. a can be found by realizing that:

b can be obtained since the hypotenuse of each right triangle is equal to half the sample distance. A little trigonometry takes care of finding b:

Eccentricity e is simply defined as:

We then plug all this into Equation (6) and get a vector of radii for all phi between pi and 2*pi (the negative half of the unit circle).

This vector is mapped to the image grid via the matlab command pol2cart and the shifted by c to bring the ellipse into line with the focii (source and receiver).

Finally, the value of the time sample is then added to these grid coordinates. The process is repeated for each source-receiver time vector.

Time Comparison

Both methods require treating each of the 4096 vectors separately and traversing the vectors one by one. Therefore, we can talk about efficiency in terms of the time required to process one vector and extrapolate this to total processing time. Without any optimizations, both methods take several seconds to process one vector, which means that total process time is on the order of 5-10 hours. So how can we reduce these times? Well, one obvious answer would be to use a different computing simulator or work with a language more optimized for this type of processing. However, for our purposes, let’s only consider optimizations to the actual algorithms.

For method one, the main factor we can control is the size of the grid. We may greatly reduce our processing time by simply searching a smaller grid. The idea is to traverse the entire grid for the first few vectors and then to only traverse the areas where the magnitude is greater than some threshold. It is possible, in this way, to decrease the area by a factor of 2-10 and thus improve the speed.

For method two, we set a threshold value and only consider time samples above this value. This is helpful in eliminating all of the zero values in the early part of the time vectors that occur before any pulse returns. We must be careful however in this optimization because if we set the threshold too high, we will be hurting our resolution and/or creating noise in the graph.

Resolution Comparison

Here are two pictures of the same image reconstructed using each method

Figure 2

(a) Image from Method 1 (b) Image from Method 2

As you can see, Method 1 shows details much more clearly because it uses a linear interpolation formula, while Method 2 rounds off the ellipse coordinates to fit them to the grid. This means that method one should be more accurate and should show less error from the discrete nature of the samples. Method 2 is faster, but method 1 shows higher resolution.

You can download copies of the matlab code for method 1, the optimizing wrapper for method 1 and method 2.