Inconsistent results from extracting rotated 1d profiles from 2d numpy array

Question

I have a 2D numpy array that represents an image (see below).

The elliptical object in the image is rotated an angle theta from the vertical y-axis. In my code, I'd like to measure the FWHM of the object at various positions (including through the one that will be maximized, the semi-major axis at angle theta).

To do this, I've used techniques from both of these questions (to extract the line at a given angle, and then use UnivariateSpline to compute the FWHM given a 1d profile): How to extract an arbitrary line of values from a numpy array? and Finding the full width half maximum of a peak

However, I've noticed that my results are inconsistent. If I extract a profile at a given angle from the entire image (which is a (641, 641) array), I get a different result for the FWHM than if I get the profile at the same angle out of a (100,100) subimage where the object is centered. I realize that since the method involves interpolation, I will likely get different values along the profiles. But what I need is a reliable and consistent way to compute these FWHMs, and the angle at which it is maximized (since right now it's giving me a different angle for subimage vs whole image). Here's my code:

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

Now, to find the angle from the y-axis at which the FWHM is maximized:

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)

Outputs: Maximized FWHM and associated position angle: 20.899 14.5

From the information on the website I linked, the pixel scale of the image is 0.275 arcseconds. So multiplying by that, the FWHM should at its maximum at a position angle of 14.5 degrees, with a value of about 5.75 arcseconds. However, Table 1 on the website clearly states that the maximized position angle from the y-axis is only 6 degrees, and that the FWHM of the major axis there is 7.36 arcseconds. So something must be wrong here.

And if I run this code again but on the whole image instead of a subimage, I get a completely different result for the angle and the FWHM. Does anyone know how I can find a more consistent (and accurate) method? Thank you!

Answer 1

Not 100% sure, but you seem to snap your line coordinates to the pixel grid which feels pretty inaccurate to me. Also, the "inverted letterbox" (white bars) might confuse your program? Maybe it would be safer to use proper 2D interpolation and clip the image, like so:

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