Vectorizing python code to numpy

Question

I have the following code snippet (for Hough circle transform):

for r in range(1, 11):
    for t in range(0, 360):
        trad = np.deg2rad(t)

        b = x - r * np.cos(trad)
        a = y - r * np.sin(trad)

        b = np.floor(b).astype('int')
        a = np.floor(a).astype('int')

        A[a, b, r-1] += 1

Where A is a 3D array of shape (height, width, 10), and height and width represent the size of a given image. My goal is to convert the snippet exclusively to numpy code.

My attempt is this:

arr_r = np.arange(1, 11)
arr_t = np.deg2rad(np.arange(0, 360))

arr_cos_t = np.cos(arr_t)
arr_sin_t = np.sin(arr_t)

arr_rcos = arr_r[..., np.newaxis] * arr_cos_t[np.newaxis, ...]
arr_rsin = arr_r[..., np.newaxis] * arr_sin_t[np.newaxis, ...]

arr_a = (y - arr_rsin).flatten().astype('int')
arr_b = (x - arr_rcos).flatten().astype('int')

Where x and y are two scalar values.

I am having trouble at converting the increment part: A[a,b,r] += 1. I thought of this: A[a,b,r] counts the number of occurrences of the pair (a,b,r), so a clue was to use a Cartesian product (but the arrays are too large).

Any tips or tricks I can use?

Thank you very much!

Edit: after filling A, I need (a,b,r) as argmax(A). The tuple (a,b,r) identifies a circle and its value in A represents the confidence value. So I want that tuple with the highest value in A. This is part of the voting algorithm from Hough circle transform: find circle parameter with unknown radius.

Is there a reason you have a dimension for r in A? It seems redundant. — Paul Panzer
– Paul Panzer, Commented Dec 28, 2017 at 20:09
@PaulPanzer - after filling A, I need (a,b,r) as argmax(A). The tuple (a,b,r) identifies a circle and its value in A represents the confidence value. So I want that tuple with the highest value in A. — Alexandru Dinu
– Alexandru Dinu, Commented Dec 28, 2017 at 20:15

Divakar · Accepted Answer · 2017-12-28 21:30:14Z

3

Method #1

Here's one way leveraging broadcasting to get the counts and update A (this assumes the a and b values computed in the intermediate steps are positive ones) -

d0,d1,d2 = A.shape    
arr_r = np.arange(1, 11)
arr_t = np.deg2rad(np.arange(0, 360))

arr_b = np.floor(x - arr_r[:,None] * np.cos(arr_t)).astype('int')
arr_a = np.floor(y - arr_r[:,None] * np.sin(arr_t)).astype('int')

idx = (arr_a*d1*d2) + (arr_b * d2) + (arr_r-1)[:,None]

A.flat[:idx.max()+1] += np.bincount(idx.ravel())
# OR A.flat += np.bincount(idx.ravel(), minlength=A.size)

Method #2

Alternatively, we could avoid bincount to replace the last step in approach #1, like so -

idx.ravel().sort()
idx.shape = (-1)
grp_idx = np.flatnonzero(np.concatenate(([True], idx[1:]!=idx[:-1],[True])))
A.flat[idx[grp_idx[:-1]]] += np.diff(grp_idx)

Improvement with numexpr

We could also leverage numexpr module for faster sine, cosine computations, like so -

import numexpr as ne

arr_r2D = arr_r[:,None]
arr_b = ne.evaluate('floor(x - arr_r2D * cos(arr_t))').astype(int)
arr_a = ne.evaluate('floor(y - arr_r2D * sin(arr_t))').astype(int)

edited Dec 28, 2017 at 21:30

answered Dec 28, 2017 at 20:27

Divakar

222k19 gold badges273 silver badges374 bronze badges

Sign up to request clarification or add additional context in comments.

8 Comments

Alexandru Dinu Over a year ago

I have to mention that those two for loops are nested under another two for loops (for x and for y). I have to compute this for every (x,y).

Divakar Over a year ago

@Alex But the question states - Where x and y are two scalar values.?

Alexandru Dinu Over a year ago

They are scalars (indices). This numpy code will be run for each (x,y).

Divakar Over a year ago

@Alex Cleaned up code. You might want to try out method#2 as well. And numexpr is worth a go as well I think.

Divakar Over a year ago

@Alex Well idx is array of flattened indices to be incremented in A and they appear as many times, the incrementing is to be done. So, we flatten idx, sort it. Then, get the start of each group of idenntical indices, that's grp_idx. Then, get count of each group of identical indices with np.diff. Finally, we use start of each group to select that flattened index corresponding to each group, index into flattened A and assign the counts.

|

Thomas · Accepted Answer · 2017-12-28 22:36:48Z

2

np.add(np.array ([arr_a, arr_b, 10]), 1)

edited Dec 28, 2017 at 22:36

answered Dec 28, 2017 at 20:09

Thomas

1,0669 silver badges19 bronze badges

3 Comments

Divakar Over a year ago

What exactly are you doing here? I am seeing typos all over and still can't see how simply indexing with those indexing arrays would work out.

Thomas Over a year ago

Sorry for the typo. Corrected it. I understood Alex wants to increment the three values in his array using numpy which the above line does. He also wrote r is always 10 so I hard-coded the value.

Divakar Over a year ago

And how are you updating A?

Collectives™ on Stack Overflow

Vectorizing python code to numpy

2 Answers 2

8 Comments

3 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

8 Comments

3 Comments

Your Answer

Sign up or log in

Post as a guest

Related