Trilateration in python

Question

I am working on a positioning system. The input I have is a dict which will give us circles of radius d1 from point(x1,y1) and so on. The output I want is an array(similar to a 2D coordinate system) in which the intersecting area is marked 1 and rest is 0. I tried this:

xsize=3000
ysize=2000
lis={(x1,y1):d1,(x2,y2):d2,(x3,y3):d3}
array=np.zeros((xsize,ysize))
for i in range(xsize-1):
    for j in range(ysize-1):
        for element in lis:
            if distance((i,j),element)<=(lis[element]):
                array[i][j]=1
            else:
                array[i][j]=0
                break
def distance(p1,p2):
    return math.sqrt((p1[0]-p2[0])**2+(p1[1]-p2[1])**2)

The only problem is that the array is large and takes way too long(no. of loops is in 10 millions), especially on a raspberry pi, otherwise this works. Is there any way to do it using openCV and an image and then draw circles to get the intersecting area faster?

It has to be python 2.x.

Please provide a [MCVE]. Currently I can't even figure out what's working that's working. — timgeb, Aug 29 '18 at 08:25
I edited it a little, did you get my problem now?(nothing wrong with the syntax just need a faster way to get the result without so many loops) @timgeb — Vedant Kandoi, Aug 29 '18 at 08:48
I guess you could create an empty black image 3000x2000 in memory in a `numpy` array. Then, for each circle in your list, create a square 2*d by 2*d black canvas (again as a `numpy` array) and draw a circle in it with `fill colour=1`. Add that square to your main image using slicing. Then at the end, if your list had 7 circles in it, create a mask of all pixels whose value is 7. Watch out for anti-aliasing, by the way. — Mark Setchell, Aug 29 '18 at 09:03
By the way, if you are testing if a distance is less than 10, it is quicker to not take the `sqrt()` (which is slow) and test if the distance squared is less than 100. — Mark Setchell, Aug 29 '18 at 09:14

w-m · Answer 1 · 2018-08-29T09:53:28.967

As you are already using numpy, try to rewrite your operations in a vectorized fashion, instead of using loops.

# choose appropriate dtype for better perf
dtype = np.float32

# take all indices in an array
indices = np.indices((ysize, xsize), dtype=dtype).T
points = np.array(list(lis.keys()), dtype=dtype)

# squared distance from all indices to all points
dist = (indices[..., np.newaxis] - points.T) ** 2
dist = dist.sum(axis=-2)

# squared circle radii
dist_thresh = np.array(list(lis.values()), dtype=dtype) ** 2

intersect = np.all(dist <= dist_thresh, axis=-1)

That's around 60x faster on my machine than the for loop version.

It is still a brute-force version, doing possibly many needless computations for all coordinates. The circles are not given in the question, so it's hard to reason about them. If they cover a relatively small area, the problem will be solved much faster (still computationally, not analytically), if a smaller area is considered. For example instead of testing all coordinates, the intersection of the bounding boxes of the circles could be used, which may reduce the computational load considerably.

score 0 · Answer 2 · edited Jun 13 '22 at 11:25

Thanks for the answers!

I also found this:

pos=np.ones((xsize,ysize))
xx,yy=np.mgrid[:xsize,:ysize]
for element in lis:
    circle=(xx-element[0])**2+(yy-element[1])**2
    pos=np.logical_and(pos,(circle<(lis[element]**2)))
    #pos&circle<(lis[element]**2 doesn't work(I read somewhere it does)

I needed this array for marking when I reached my destination or not.

if pos[dest[0]][dest[1]]==1 #Reached

Trilateration in python

2 Answers2