CODEWITHSUNDEEP
pythonnumpymultidimensional-array
Philipp
Philipp

Reputation: 2516

Displaying a slice from a 3d ndarray

I have a Hough Transform that tries to find circles where both center and radius is unknown so the accumulator space is 3 dimensional.

DIMENSION = 200
R_DIMENSION = 200

xbins = np.linspace(-0.5,0.5,DIMENSION)
ybins = np.linspace(-0.5,0.5,DIMENSION)
rbins = np.linspace(0,0.5, R_DIMENSION)

x,y,r = np.broadcast_arrays( xbins[...,np.newaxis,np.newaxis], \
                               ybins[np.newaxis,...,np.newaxis], \
                               rbins[np.newaxis,np.newaxis,...])
score = 0
circle_counter = 0
while True:
  weights = np.zeros( (DIMENSION, DIMENSION, R_DIMENSION))
  for x0,y0 in data['allPoints']:
    s = 0.001
    eta = (x-x0)**2 + (y-y0)**2 - r**2
    weights += 1. / ( sqrt( 2 * sconst.pi ) * s ) * np.exp( -( eta ** 2 )\
                                                              / ( 2 * s ** 2 ) )
  index = np.argmax( weights )
  ii,jj,rr = np.unravel_index( index, (DIMENSION, DIMENSION, R_DIMENSION))
  score = weights[ii][jj][rr]
  if score < 200:
    break

Now if I want to visualise the x,y space for rr for some reason I don't get the plane with the highest score BUT using the unravelled indices for determining the actual radius and center gives me the right result when I plot it.

circle = {}
circle['center'] = (xbins[ii], ybins[jj])
circle['radius'] = rbins[rr]
circles.append(circle)

plt.imshow(weights[:][:][rr])
plt.colorbar()
plt.show()

So my question is do I misunderstand how I display the x,y plane for a given radius index?

Here are two pictures of the visualisation. The 2nd one contains the actual maximum. What confuses me in both pictures is that the lines are warped meaning it doesn't seem like it is for a fixed radius.

first picture was created with imshow(weights[:][:][rr]) and the second one is a from a series

for r_i in range(R_DIMENSION):
  imshow(weights[:][:][r_i])

and then I find the one that contains the highest score. Slice with 3rd dimension fixed with <code>rr</code>

real maximum

and actually I expect something like this (this is from a 2D hough transform where I know the radius): enter image description here

Upvotes: 1

Views: 133

Answers (1)

Philipp
Philipp

Reputation: 2516

OK, I found a solution, but if someone could explain it to me that would be great.

What I had to change is how I broadcast the arrays:

x,y,r = np.broadcast_arrays( xbins[np.newaxis,...,np.newaxis], \
                             ybins[np.newaxis,np.newaxis,...], \
                             rbins[...,np.newaxis,np.newaxis])

and of course then changing the order of the indices/dimensions where I use them. 

while True:
  weights = np.zeros( (R_DIMENSION, DIMENSION, DIMENSION))
  for x0,y0 in data['allPoints']:
    s = 0.001
    eta = (x-x0)**2 + (y-y0)**2 - r**2
    weights += 1. / ( sqrt( 2 * sconst.pi ) * s ) * np.exp( -( eta ** 2 )\
                                                            / ( 2 * s ** 2 ) )
  index = np.argmax( weights )
  rr,ii,jj = np.unravel_index( index, (R_DIMENSION, DIMENSION, DIMENSION))
  score = weights[rr][ii][jj]
  if score < 200:
    print 'finished after %s circle(s) found' % circle_counter
    break

The animation of different radiuses for a x,y plane then looks like this.

Upvotes: 1

Related Questions