0

My goal is to compute the FWHM of a star profile.

I've an image and a single star in it, as input. Each pixel, at position (x, y), has an intensity associated between 0 and 1.

My idea was to compute the standard deviation of the entire dataset and then use the following formula:

f(x, y)=[1/(√(2π)σ)]exp(-[(x - X)^2 + (y - Y)^2]/2σ^2])

solving the equation:

fmax/2=1/[2√(2π)σ]=[1/(√(2π)σ)]exp(-[(x - X)^2 + (y - Y)^2]/2σ^2]) =>

FWHM=2σ√(2ln2)

With this approach I'm not getting the expected result looking at the data.

Is there something I'm missing? Any other suggestion?

Minal Chauhan
  • 6,025
  • 8
  • 21
  • 41
rcolombari
  • 79
  • 10

1 Answers1

0

You might get better results using interpolation, so that you can achieve subpixel precision. And since your data is between 0-1, with your centroid having an intensity of 1, you could extract the profile(s) through that point. First, import your image as a 2D numpy array if it isn't in that form already. If you're starting from a FITS file:

from astropy.io import fits
filename = 'your_image_file.fits'
image = fits.getdata(filename)

To extract your profile and get the FWHM:

import numpy as np
from scipy.interpolate import UnivariateSpline

def profiles(image):
    ypix, xpix = np.where(image==1)
    x = np.take(image, ypix[0], axis=0)
    y = np.take(image, xpix[0], axis=1)

    return x, y #these are the horizontal and vertical profiles through the star's centroid

def interpolate_width(axis):
    half_max = 1/2
    x = np.linspace(0, len(axis), len(axis))

    # Do the interpolation
    spline = UnivariateSpline(x, axis-half_max, s=0)
    r1, r2 = spline.roots()

    return r2-r1 #this is the FWHM along the specified axis

horizontal, vertical = profiles(image)
fwhm_x = interpolate_width(horizontal)
fwhm_y = interpolate_width(vertical)

This is assuming that the star is not rotated -- or if it is, that you just want the FWHMs along the horizontal and vertical axes. If the star is rotated with respect to the horizontal, and you want the FWHMs along the semi-major and semi-minor axes, you'll have to extract the profile by taking the line segment connected by two points. Then you can get the FWHM the same way with the interpolate_width function. See here for the profile extraction portion of this approach: How to extract an arbitrary line of values from a numpy array?

curious_cosmo
  • 1,184
  • 1
  • 18
  • 36