Multiliteration implementation with inaccurate distance data

Question

I am trying to create an android smartphone application which uses Apples iBeacon technology to determine the current indoor location of itself. I already managed to get all available beacons and calculate the distance to them via the rssi signal.

Currently I face the problem, that I am not able to find any library or implementation of an algorithm, which calculates the estimated location in 2D by using 3 (or more) distances of fixed points with the condition, that these distances are not accurate (which means, that the three "trilateration-circles" do not intersect in one point).

I would be deeply grateful if anybody can post me a link or an implementation of that in any common programming language (Java, C++, Python, PHP, Javascript or whatever). I already read a lot on stackoverflow about that topic, but could not find any answer I were able to convert in code (only some mathematical approaches with matrices and inverting them, calculating with vectors or stuff like that).

EDIT

I thought about an own approach, which works quite well for me, but is not that efficient and scientific. I iterate over every meter (or like in my example 0.1 meter) of the location grid and calculate the possibility of that location to be the actual position of the handset by comparing the distance of that location to all beacons and the distance I calculate with the received rssi signal.

Code example:

public Location trilaterate(ArrayList<Beacon> beacons, double maxX, double maxY)
{
    for (double x = 0; x <= maxX; x += .1)
    {
        for (double y = 0; y <= maxY; y += .1)
        {
            double currentLocationProbability = 0;
            for (Beacon beacon : beacons)
            {
                // distance difference between calculated distance to beacon transmitter
                // (rssi-calculated distance) and current location:
                // |sqrt(dX^2 + dY^2) - distanceToTransmitter|
                double distanceDifference = Math
                    .abs(Math.sqrt(Math.pow(beacon.getLocation().x - x, 2)
                                   + Math.pow(beacon.getLocation().y - y, 2))
                         - beacon.getCurrentDistanceToTransmitter());
                // weight the distance difference with the beacon calculated rssi-distance. The
                // smaller the calculated rssi-distance is, the more the distance difference
                // will be weighted (it is assumed, that nearer beacons measure the distance
                // more accurate)
                distanceDifference /= Math.pow(beacon.getCurrentDistanceToTransmitter(), 0.9);
                // sum up all weighted distance differences for every beacon in
                // "currentLocationProbability"
                currentLocationProbability += distanceDifference;
            }
            addToLocationMap(currentLocationProbability, x, y);
            // the previous line is my approach, I create a Set of Locations with the 5 most probable locations in it to estimate the accuracy of the measurement afterwards. If that is not necessary, a simple variable assignment for the most probable location would do the job also
        }
    }
    Location bestLocation = getLocationSet().first().location;
    bestLocation.accuracy = calculateLocationAccuracy();
    Log.w("TRILATERATION", "Location " + bestLocation + " best with accuracy "
                           + bestLocation.accuracy);
    return bestLocation;
}

Of course, the downside of that is, that I have on a 300m² floor 30.000 locations I had to iterate over and measure the distance to every single beacon I got a signal from (if that would be 5, I do 150.000 calculations only for determine a single location). That's a lot - so I will let the question open and hope for some further solutions or a good improvement of this existing solution in order to make it more efficient.

Of course it has not to be a Trilateration approach, like the original title of this question was, it is also good to have an algorithm which includes more than three beacons for the location determination (Multilateration).

Hagbut

Answer 1

If the current approach is fine except for being too slow, then you could speed it up by recursively subdividing the plane. This works sort of like finding nearest neighbors in a kd-tree. Suppose that we are given an axis-aligned box and wish to find the approximate best solution in the box. If the box is small enough, then return the center.

Otherwise, divide the box in half, either by x or by y depending on which side is longer. For both halves, compute a bound on the solution quality as follows. Since the objective function is additive, sum lower bounds for each beacon. The lower bound for a beacon is the distance of the circle to the box, times the scaling factor. Recursively find the best solution in the child with the lower lower bound. Examine the other child only if the best solution in the first child is worse than the other child's lower bound.

Most of the implementation work here is the box-to-circle distance computation. Since the box is axis-aligned, we can use interval arithmetic to determine the precise range of distances from box points to the circle center.

_{P.S.: Math.hypot is a nice function for computing 2D Euclidean distances.}

Anglocatholic

Answer 2

Instead of taking confidence levels of individual beacons into account, I would instead try to assign an overall confidence level for your result after you make the best guess you can with the available data. I don't think the only available metric (perceived power) is a good indication of accuracy. With poor geometry or a misbehaving beacon, you could be trusting poor data highly. It might make better sense to come up with an overall confidence level based on how well the perceived distance to the beacons line up with the calculated point assuming you trust all beacons equally.

I wrote some Python below that comes up with a best guess based on the provided data in the 3-beacon case by calculating the two points of intersection of circles for the first two beacons and then choosing the point that best matches the third. It's meant to get started on the problem and is not a final solution. If beacons don't intersect, it slightly increases the radius of each up until they do meet or a threshold is met. Likewise, it makes sure the third beacon agrees within a settable threshold. For n-beacons, I would pick 3 or 4 of the strongest signals and use those. There are tons of optimizations that could be done and I think this is a trial-by-fire problem due to the unwieldy nature of beaconing.

import math

beacons = [[0.0,0.0,7.0],[0.0,10.0,7.0],[10.0,5.0,16.0]] # x, y, radius

def point_dist(x1,y1,x2,y2):
    x = x2-x1
    y = y2-y1
    return math.sqrt((x*x)+(y*y))

# determines two points of intersection for two circles [x,y,radius]
# returns None if the circles do not intersect
def circle_intersection(beacon1,beacon2):
    r1 = beacon1[2]
    r2 = beacon2[2]
    dist = point_dist(beacon1[0],beacon1[1],beacon2[0],beacon2[1])
    heron_root = (dist+r1+r2)*(-dist+r1+r2)*(dist-r1+r2)*(dist+r1-r2)
    if ( heron_root > 0 ):
        heron = 0.25*math.sqrt(heron_root)
        xbase = (0.5)*(beacon1[0]+beacon2[0]) + (0.5)*(beacon2[0]-beacon1[0])*(r1*r1-r2*r2)/(dist*dist)
        xdiff = 2*(beacon2[1]-beacon1[1])*heron/(dist*dist) 
        ybase = (0.5)*(beacon1[1]+beacon2[1]) + (0.5)*(beacon2[1]-beacon1[1])*(r1*r1-r2*r2)/(dist*dist)
        ydiff = 2*(beacon2[0]-beacon1[0])*heron/(dist*dist) 
        return (xbase+xdiff,ybase-ydiff),(xbase-xdiff,ybase+ydiff)
    else:
        # no intersection, need to pseudo-increase beacon power and try again
        return None

# find the two points of intersection between beacon0 and beacon1
# will use beacon2 to determine the better of the two points
failing = True
power_increases = 0
while failing and power_increases < 10:
    res = circle_intersection(beacons[0],beacons[1])
    if ( res ):
        intersection = res
    else:
        beacons[0][2] *= 1.001
        beacons[1][2] *= 1.001
        power_increases += 1
        continue
    failing = False

# make sure the best fit is within x% (10% of the total distance from the 3rd beacon in this case)
# otherwise the results are too far off
THRESHOLD = 0.1

if failing:
    print 'Bad Beacon Data (Beacon0 & Beacon1 don\'t intersection after many "power increases")'
else:
    # finding best point between beacon1 and beacon2
    dist1 = point_dist(beacons[2][0],beacons[2][1],intersection[0][0],intersection[0][1])
    dist2 = point_dist(beacons[2][0],beacons[2][1],intersection[1][0],intersection[1][1])
    if ( math.fabs(dist1-beacons[2][2]) < math.fabs(dist2-beacons[2][2]) ):
        best_point = intersection[0]
        best_dist = dist1
    else:
        best_point = intersection[1]
        best_dist = dist2
    best_dist_diff = math.fabs(best_dist-beacons[2][2])
    if best_dist_diff < THRESHOLD*best_dist:
        print best_point
    else:
        print 'Bad Beacon Data (Beacon2 distance to best point not within threshold)'

If you want to trust closer beacons more, you may want to calculate the intersection points between the two closest beacons and then use the farther beacon to tie-break. Keep in mind that almost anything you do with "confidence levels" for the individual measurements will be a hack at best. Since you will always be working with very bad data, you will defintiely need to loosen up the power_increases limit and threshold percentage.

Formicary

Answer 3

You have 3 points : A(xA,yA,zA), B(xB,yB,zB) and C(xC,yC,zC), which respectively are approximately at dA, dB and dC from you goal point G(xG,yG,zG). Let's say cA, cB and cC are the confidence rate ( 0 < cX <= 1 ) of each point. Basically, you might take something really close to 1, like {0.95,0.97,0.99}. If you don't know, try different coefficient depending of distance avg. If distance is really big, you're likely to be not very confident about it.

Here is the way i'll do it :

var sum = (cA*dA) + (cB*dB) + (cC*dC);
dA = cA*dA/sum;
dB = cB*dB/sum;
dC = cC*dC/sum;

xG = (xA*dA) + (xB*dB) + (xC*dC);
yG = (yA*dA) + (yB*dB) + (yC*dC);
xG = (zA*dA) + (zB*dB) + (zC*dC);

Basic, and not really smart but will do the job for some simple tasks.

EDIT

You can take any confidence coef you want in [0,inf[, but IMHO, restraining at [0,1] is a good idea to keep a realistic result.

Selwyn

I take me a while to understand your approach. I think the problem is, that it tries to find a point within the 3 given points. For example, if I have the following setup of points, ranges (all with confidence 1.0, in this example measurement is exact): [{x: 0, y: 0, z: 0, d: 3}, {x: 1, y: 3, z: 0, d: 1}, {x: 4, y: 0, z: 0, d: 5}] and the goal point is on {x: 0, y: 3, z: 0}. Unfortunately your algorithm does not give any point close to that: {x: 2.3333333333333335, y: 0.3333333333333333, z: 0} is the result of the calculation. Do you have any idea for an improvement or did I miss something? – Hagbut