Search code examples
pythonnumpymaskcalculus

Compute sum of image pixels along a circular path


I have an image where I am trying to compute the line integral (sum) along a circular path. My idea to approach this is:

  1. Compute the circular path to sum over
  2. Mask the image based on the path, where everything except whatever pixels coincide with the path are zeroed out.
  3. Sum over all image pixel

I am currently stuck between steps one and two, where I can't figure out how to generate a circle on the same grid as the image.

In code:

from scipy.stats import multivariate_normal

radius = 2

# Draw arbitrary image 
x, y = np.mgrid[-5:5:.1, -5:5:.1]
img = multivariate_normal.pdf(np.dstack((x, y)), cov=[[1, 0.7], [0.7, 1]])

# generate circle with desired radius 
circle = radius*np.exp(1j*np.linspace(-np.pi, np.pi, 100))

pyplot.pcolormesh(x, y, img)
pyplot.plot(np.real(circle), np.imag(circle), '-w')
pyplot.show()

Which gives: enter image description here

Question:

How to use the circle to mask the image pixels coinciding with this circle?


Solution

  • Here is an alternative way of calculating the integral: It uses interpolation so the image becomes a function defined on a rectangle and then computes the path integral using a standard integral solver.

    from scipy.integrate import quad
    from scipy.interpolate import RectBivariateSpline
    from scipy.stats import multivariate_normal
    import numpy as np
    
    x, y = np.ogrid[-5:5:.1, -5:5:.1]
    img = multivariate_normal.pdf(np.dstack(np.broadcast_arrays(x, y)),
                                  cov=[[1, 0.7], [0.7, 1]])
    
    f = RectBivariateSpline(x.ravel(), y.ravel(), img)
    
    radius, centerx, centery = 3.0, 1.0, -1.5
    def integrand(rad):
        return f(centerx+radius*np.cos(rad), centery+radius*np.sin(rad))
    
    def true_integrand(rad):
        return multivariate_normal(cov=[[1, 0.7], [0.7, 1]]).pdf(
            (centerx+radius*np.cos(rad), centery+radius*np.sin(rad)))
    
    print(quad(integrand, -np.pi, np.pi))
    print(quad(true_integrand, -np.pi, np.pi))
    

    Output:

    (0.07985467350026378, 1.3411796499850778e-08)
    (0.07985453947958436, 4.006916325573184e-11)