Вычислить сумму пикселей изображения по круговой траектории

У меня есть изображение, где я пытаюсь вычислить линейный интеграл (сумма) по круговой траектории. Моя идея подойти к этому:

  1. Вычислить круговой путь для суммирования по
  2. Маскируйте изображение на основе пути, где все, кроме пикселей, совпадают с путем, обнуляются.
  3. Сумма по всему пикселю изображения

В настоящее время я застрял между шагами один и два, где я не могу понять, как создать круг на той же сетке, что и изображение.

В коде:

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

Который дает:

Вопрос:

Как использовать круг для маскировки пикселей изображения, совпадающих с этим кругом?

1 ответ

Решение

Вот альтернативный способ вычисления интеграла: он использует интерполяцию, поэтому изображение становится функцией, определенной в прямоугольнике, а затем вычисляет интеграл по траектории с использованием стандартного решателя интегралов.

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

Выход:

(0.07985467350026378, 1.3411796499850778e-08)
(0.07985453947958436, 4.006916325573184e-11)
Другие вопросы по тегам