屏蔽位于两条曲线之间的图像（np.ndarray）部分

Question

从astropy借用这个例子：

import numpy as np
from matplotlib import pyplot as plt
from astropy.io import fits
from astropy.wcs import WCS
from astropy.utils.data import download_file

fits_file = 'http://data.astropy.org/tutorials/FITS-images/HorseHead.fits'
image_file = download_file(fits_file, cache=True)
hdu = fits.open(image_file)[0]
wcs = WCS(hdu.header)

fig = plt.figure()
ax = fig.add_subplot(111)
plt.imshow(hdu.data, origin='lower', cmap='cubehelix')
plt.xlabel('X')
plt.ylabel('Y')

x_array = np.arange(0, 1000)
line_1 = 1 * x_array + 20 * np.sin(0.05*x_array)
line_2 = x_array - 100 + 20 * np.sin(0.05*x_array)

plt.plot(x_array, line_1, color='red')
plt.plot(x_array, line_2, color='red')

ax.set_xlim(0, hdu.shape[1])
ax.set_ylim(0, hdu.shape[0])

plt.show()

我想计算位于两条曲线之间的像素中值（例如在 y 方向）：

我相信聪明的做法是为感兴趣的区域创建一个蒙版。

有没有办法在不循环图像像素的情况下生成此蒙版？

编辑 1：修改后的问题提高了理解

编辑 2：我更改了示例以更好地表示问题标题（曲线而不是直线）

Answer 1

生成这样一个掩码的一种相当直接的方法是使用numpy.mgrid thingy，它本质上为您提供了 x 和 y 坐标数组（在那种 2D 情况下），然后可以使用该等式计算掩码的行。

编辑：只要您可以使用方程式（例如f(x,y)<0，其中 f 是您想要的任何函数）或这些方程式的组合来表达您的掩码，您就可以做任何您想做的事情。以下是您的新面具以及一些额外艺术作品的示例：

import numpy as np
import matplotlib.pyplot as plt

y,x=np.mgrid[0:1000,0:1000]

#a wiggly ramp
plt.subplot(221)
mask0=(y < x + 20 * np.sin(0.05*x)) & (y > x - 100 + 20 * np.sin(0.05*x))
plt.imshow(mask0,origin='lower',cmap='gray')

#a ramp
plt.subplot(222)
mask1=(y < x) & (x < y+100)
plt.imshow(mask1,origin='lower',cmap='gray')

#a disk
plt.subplot(223)
mask2=(200**2>(x-500)**2+(y-500)**2)
plt.imshow(mask2,origin='lower',cmap='gray')

#a ying-yang attempt
plt.subplot(224)
mask3= (mask2 & (0 < np.sin(3.14*x/250)*100 + 500 - y) & (30**2 < (x-620)**2+(y-500)**2) )| (30**2 > (x-380)**2+(y-500)**2)
plt.imshow(mask3,origin='lower',cmap='gray')

plt.show()

输出：

Answer 2

下面的代码 sn-p 创建一个掩码，将掩码外的所有值设置为 nan，然后使用 NumPy 的nanmedian 沿着指示的方向计算所需的数量。

import numpy as np

import matplotlib.pyplot as plt
from matplotlib import gridspec

from skimage.measure import grid_points_in_poly


# Create test image
N = 900
image = np.sin(np.linspace(0, 2 * np.pi, N * N).reshape((N, N)) ** 2)

# Define the mask polygon
poly = [[N, 0],
        [0, N],
        [100, N],
        [N, 100]]

# Create the mask (True is inside the polygon)
mask = grid_points_in_poly(image.shape, poly)

# Set everything outside the mask to nan
masked_image = image.copy()
masked_image[~mask] = np.nan

# Perform the required operation
row_med = np.nanmedian(masked_image, axis=1)

# The rest of the code is to visualize the result
fig = plt.figure(figsize=(15, 10))
gs = gridspec.GridSpec(1, 3, width_ratios=(1, 1, 1/8))
ax0 = plt.subplot(gs[0])
ax1 = plt.subplot(gs[1])
ax2 = plt.subplot(gs[2])

ax0.imshow(image, cmap='gray')
ax0.set_title('Input image')
ax0.set_xlim(0, N)

ax1.imshow(masked_image, cmap='gray')
ax1.set_title('Masked image')
ax1.set_xlim(0, N)

ax2.plot(row_med, np.arange(N))
ax2.set_ylim([N, 0])
ax2.set_xlim([-1.5, 1.5])
ax2.set_title('Row median')

plt.tight_layout()
plt.show()