3

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?

Gerges
  • 6,269
  • 2
  • 22
  • 44
  • You can generate circle coordinates using Midpoint circle alg. [link](https://en.wikipedia.org/wiki/Midpoint_circle_algorithm) and then: `np.sum(image[circle_coordinates])` – Brenlla Apr 17 '18 at 17:51
  • @Brenlla I think to be correct one has to account for path lengths which are different for different pixels. For example the pixel at `(1/sqrt(2), 1/sqrt(2))` should be weighted `sqrt(2)` times more than the pixel at `(1, 0)` (assuming `radius=1`). – Paul Panzer Apr 17 '18 at 18:42

1 Answers1

3

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)
Paul Panzer
  • 51,835
  • 3
  • 54
  • 99
  • 1
    Thank you, this is a clever way to go around the computation! – Gerges Apr 17 '18 at 21:23
  • I tried testing this by using `img = np.ones((max(x.shape), max(y.shape)))`, so I would compute the circumference, which should be 2 * pi * radius, but this only gives me 2 * pi. Any ideas? I can post this as a new question if there is no quick answer. – Gerges Apr 17 '18 at 21:45
  • Never mind, I interpreted this wrongly, it does give the expected answer, sorry for spamming your answer :). – Gerges Apr 17 '18 at 22:36