从2d numpy数组中提取旋转一维轮廓的结果不一致

Question

我有一个表示图像的2D numpy数组（见下文）。

图像中的椭圆物体从垂直y轴旋转theta。在我的代码中，我想测量物体在不同位置的半高宽（包括通过一个将被最大化的位置，角度为theta的半长轴）。你知道吗

为了做到这一点，我使用了这两个问题的技巧（以给定的角度提取线，然后使用单变量样条线来计算给定一维剖面的半高宽）：How to extract an arbitrary line of values from a numpy array?和Finding the full width half maximum of a peak

然而，我注意到我的结果不一致。如果我从整个图像（这是一个（641641）数组）中提取一个给定角度的轮廓，我得到的半高宽结果与从一个（100100）子图像中得到相同角度的轮廓的结果不同。我意识到，由于该方法涉及插值，我可能会得到沿轮廓不同的值。但我需要的是一个可靠和一致的方法来计算这些FWHM，以及它最大化的角度（因为现在它给了我一个不同的角度，子图像和整个图像）。这是我的密码：

import numpy as np
from scipy.interpolate import UnivariateSpline

def profiles(image, theta):
    max_y, max_x = np.where(image==1) #the image has been normalized so that the highest intensity point=1
    max_y, max_x = max_y[0], max_x[0]

    subimage = image[max_y-50:max_y+50, max_x-50:max_x+50] # now the maximum is exactly centered at (50, 50)

    # Set up to extract profiles (do two legs to ensure it passes through centroid)
    x_maj_1, y_maj_1 = 50-50*np.tan(theta), 0
    mid_x, mid_y = 50, 50
    x_maj_2, y_maj_2 = 50+50*np.tan(theta), 99

    # Form axes
    length_maj = int(round(np.hypot(x_maj_2-x_maj_1, y_maj_2-y_maj_1)))
    x1, y1 = np.linspace(x_maj_1, mid_x, length_maj//2), np.linspace(y_maj_1, mid_y, length_maj//2)
    x2, y2 = np.linspace(mid_x, x_maj_2, length_maj//2)[1:], np.linspace(mid_y, y_maj_2, length_maj//2)[1:]

    # Concatenate legs
    x_maj = np.concatenate((x1, x2), axis=0)
    y_maj = np.concatenate((y1, y2), axis=0)

    # Get profile
    z_maj = subimage[y_maj.astype(np.int), x_maj.astype(np.int)]
    return z_maj

def interpolate_width(axis):
    half_max = 1/2
    x = np.arange(0, len(axis))
    spline = UnivariateSpline(x, axis-half_max, s=0)
    r1, r2 = spline.roots()
    return r2-r1 #FWHM in pixel units

现在，从y轴上找出半高宽最大的角度：

thetas = np.arange(0, 45, 0.5)
widths = []
for theta in thetas:
    theta = np.deg2rad(theta)
    z = profiles(image, theta)
    width = interpolate_width(z)
    widths.append(width)

fwhm_maj = max(widths)
angle_arg = np.array(widths).argmax()
angle_max = thetas[angle_arg]
print('Maximized FWHM and associated position angle:', fwhm_maj, angle_max)

输出：Maximized FWHM and associated position angle: 20.899 14.5

从我链接的网站上的信息来看，图像的像素比例是0.275弧秒。因此乘以这个，半高宽应该在14.5度的位置角达到最大值，值约为5.75弧秒。然而，网站上的表1清楚地表明，从y轴的最大位置角只有6度，长轴的半高宽为7.36弧秒。所以这里一定出了什么问题。你知道吗

如果我再次运行这个代码，但是在整个图像上，而不是在子图像上，我会得到一个完全不同的角度和半高宽的结果。有人知道我怎样才能找到一个更一致（和准确）的方法吗？谢谢您！你知道吗

Answer 1

不是100%确定，但你似乎把你的线坐标捕捉到像素网格上，这让我觉得很不准确。此外，“倒排信箱”（白色条）可能会混淆您的程序？也许使用适当的2D插值和剪辑图像会更安全，例如：

import numpy as np
from scipy import ndimage, interpolate, optimize

img = ndimage.io.imread('VJQwQ.png')
# b&w
img = img[..., 0]
# cut "white letterbox"
img = img[:, np.where(np.any(img!=255, axis=0))[0]]

# setup interpolator for off grid pixel values
x, y = img.shape
x, y = (np.arange(z) - (z-1)/2 for z in (x, y))
intr = interpolate.RectBivariateSpline(x, y, img, kx=3, ky=3)
s = np.arange(-50, 51)

# setup line sections
# for simplicity we assume the peak is in the center
def at_angle(phi):
    def f(s, shift=0):
        x, y = np.cos(phi)*s, np.sin(phi)*s
        return intr(x, y, grid=False) - shift
    return f

# example
phi = np.pi/3
f = at_angle(phi)
mx = f(0)
left = optimize.brentq(f, -50, 0, (mx/2,))
right = optimize.brentq(f, 0, 50, (mx/2,))

# plot it
from matplotlib import pylab
pylab.plot(s, f(s))
pylab.plot([left, left], [0, mx])
pylab.plot([right, right], [0, mx])
pylab.plot([-50, 50], [mx/2, mx/2])
pylab.savefig('tst.png')