Take data from a circle in python

Question

I am looking into how the intensity of a ring changes depending on angle. Here is an example of an image:

What I would like to do is take a circle of values from within the center of that doughnut and plot them vs angle. What I'm currently doing is using scipy.ndimage.interpolation.rotate and taking slices radially through the ring, and extracting the maximum of the two peaks and plotting those vs angle.

    crop = np.ones((width,width)) #this is my image
    slices = np.arange(0,width,1)
    stack = np.zeros((2*width,len(slices)))
    angles = np.linspace(0,2*np.pi,len(crop2))

    for j in range(len(slices2)): # take slices
           stack[:,j] = rotate(crop,slices[j],reshape=False)[:,width]

However I don't think this is doing what I'm actually looking for. I'm mostly struggling with how to extract the data I want. I have also tried applying a mask which looks like this;

to the image, but then I don't know how to get the values within that mask in the correct order (ie. in order of increasing angle 0 - 2pi)

Any other ideas would be of great help!

score 3 · Accepted Answer · answered Apr 14 '16 at 16:54

3

I made a different input image to help verifying correctness:

import numpy as np
import scipy as sp
import scipy.interpolate
import matplotlib.pyplot as plt

# Mock up an image.
W = 100
x = np.arange(W)
y = np.arange(W)
xx,yy = np.meshgrid(x,y)

image = xx//5*5 + yy//5*5
image = image / np.max(image)  # scale into [0,1]

plt.imshow(image, interpolation='nearest', cmap='gray')
plt.show()

To sample values from circular paths in the image, we first build an interpolator because we want to access arbitrary locations. We also vectorize it to be faster.
Then, we generate the coordinates of N points on the circle's circumference using the parametric definition of the circle x(t) = sin(t), y(t) = cos(t).
N should be at least twice the circumference (Nyquist–Shannon sampling theorem).

interp = sp.interpolate.interp2d(x, y, image)
vinterp = np.vectorize(interp)

for r in (15, 30, 45):    # radii for circles around image's center
    xcenter = len(x)/2
    ycenter = len(y)/2
    arclen = 2*np.pi*r
    angle = np.linspace(0, 2*np.pi, arclen*2, endpoint=False)
    value = vinterp(xcenter + r*np.sin(angle),
                    ycenter + r*np.cos(angle))
    plt.plot(angle, value, label='r={}'.format(r))

plt.legend()
plt.show()

answered Apr 14 '16 at 16:54

Norman

1,765
1
18
24

Hi, thanks so much for answering. The code looks great, however I get what you have here when I take out the line. image = image/np.amax(image) then it seems to work as usual! If not, I end up with the original image being entirely black apart from one small white square in the bottom right hand order. I am a little confused as to how the scipy interp2d works - I import an image using scipy.misc's imread function, so I'm not sure how to split that up into the required x y and z inputs for interp2d? – arbitraryknowledge Apr 14 '16 at 17:56
1

Ah, I wrote Python3 code, and you've got Python 2.7, right? Try `image = image * 1.0 / np.max(image)`. Also change to `xcenter = 0.5*len(x) `, etc. It's because of [integer division](http://stackoverflow.com/q/21316968/1292641). – Norman Apr 14 '16 at 18:09
1

Your image from `imread()` is an ndarray. Try `h, w = image.shape; sp.interpolate.interp2d(np.arange(w), np.arange(h), image)`. – Norman Apr 14 '16 at 18:21
That did the trick, thank you very much! I was forgetting that x and y have to be input 2d arrays. Will mark everything as solved :) – arbitraryknowledge Apr 14 '16 at 18:27
1

Correction again because of integer division: `sp.interpolate.interp2d(np.arange(w*1.0), np.arange(h*1.0), image)`. Or maybe not. I don't know whether interp2d() would handle an int array correctly. – Norman Apr 14 '16 at 18:29
1

You mean "x and y have to be input *1D* arrays". – Norman Apr 14 '16 at 18:31

Take data from a circle in python

1 Answers1