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!

Robson

Answer 1

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()

Trahern

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? – Robson

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. – Trahern

Your image from imread() is an ndarray. Try h, w = image.shape; sp.interpolate.interp2d(np.arange(w), np.arange(h), image). – Trahern

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 :) – Robson

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. – Trahern

You mean "x and y have to be input 1D arrays". – Trahern