Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def psf_norm_2d(psf, fwhm, threshold, mask_core, full_output, verbose):
""" 2d case """
# we check if the psf is centered and fix it if needed
cy, cx = frame_center(psf, verbose=False)
xcom, ycom = photutils.centroid_com(psf)
if not (np.allclose(cy, ycom, atol=1e-2) or
np.allclose(cx, xcom, atol=1e-2)):
# first we find the centroid and put it in the center of the array
centry, centrx = fit_2d(psf, full_output=False, debug=False)
shiftx, shifty = centrx - cx, centry - cy
psf = frame_shift(psf, -shifty, -shiftx, imlib=imlib,
interpolation=interpolation)
for _ in range(2):
centry, centrx = fit_2d(psf, full_output=False, debug=False)
cy, cx = frame_center(psf, verbose=False)
shiftx, shifty = centrx - cx, centry - cy
psf = frame_shift(psf, -shifty, -shiftx, imlib=imlib,
interpolation=interpolation)
# we check whether the flux is normalized and fix it if needed
imside = array.shape[0]
psf_subimage, suby, subx = get_square(array, min(cropsize, imside),
ceny, cenx, position=True)
else:
psf_subimage = array.copy()
if threshold:
_, clipmed, clipstd = sigma_clipped_stats(psf_subimage, sigma=2)
indi = np.where(psf_subimage <= clipmed + sigfactor * clipstd)
subimnoise = np.random.randn(psf_subimage.shape[0],
psf_subimage.shape[1]) * clipstd
psf_subimage[indi] = subimnoise[indi]
# Creating the 2D Gaussian model
init_amplitude = np.ptp(psf_subimage)
xcom, ycom = photutils.centroid_com(psf_subimage)
gauss = models.Gaussian2D(amplitude=init_amplitude, theta=theta,
x_mean=xcom, y_mean=ycom,
x_stddev=fwhmx * gaussian_fwhm_to_sigma,
y_stddev=fwhmy * gaussian_fwhm_to_sigma)
# Levenberg-Marquardt algorithm
fitter = fitting.LevMarLSQFitter()
y, x = np.indices(psf_subimage.shape)
fit = fitter(gauss, x, y, psf_subimage)
if crop:
mean_y = fit.y_mean.value + suby
mean_x = fit.x_mean.value + subx
else:
mean_y = fit.y_mean.value
mean_x = fit.x_mean.value
fwhm_y = fit.y_stddev.value*gaussian_sigma_to_fwhm
imside = array.shape[0]
psf_subimage, suby, subx = get_square(array, min(cropsize, imside),
ceny, cenx, position=True)
else:
psf_subimage = array.copy()
if threshold:
_, clipmed, clipstd = sigma_clipped_stats(psf_subimage, sigma=2)
indi = np.where(psf_subimage <= clipmed + sigfactor * clipstd)
subimnoise = np.random.randn(psf_subimage.shape[0],
psf_subimage.shape[1]) * clipstd
psf_subimage[indi] = subimnoise[indi]
# Creating the 2D Moffat model
init_amplitude = np.ptp(psf_subimage)
xcom, ycom = photutils.centroid_com(psf_subimage)
moffat = models.Moffat2D(amplitude=init_amplitude, x_0=xcom, y_0=ycom,
gamma=fwhm / 2., alpha=1)
# Levenberg-Marquardt algorithm
fitter = fitting.LevMarLSQFitter()
y, x = np.indices(psf_subimage.shape)
fit = fitter(moffat, x, y, psf_subimage)
if crop:
mean_y = fit.y_0.value + suby
mean_x = fit.x_0.value + subx
else:
mean_y = fit.y_0.value
mean_x = fit.x_0.value
fwhm = fit.fwhm
amplitude = fit.amplitude.value
imside = array.shape[0]
psf_subimage, suby, subx = get_square(array, min(cropsize, imside),
ceny, cenx, position=True)
else:
psf_subimage = array.copy()
if threshold:
_, clipmed, clipstd = sigma_clipped_stats(psf_subimage, sigma=2)
indi = np.where(psf_subimage <= clipmed + sigfactor * clipstd)
subimnoise = np.random.randn(psf_subimage.shape[0],
psf_subimage.shape[1]) * clipstd
psf_subimage[indi] = subimnoise[indi]
# Creating the 2d Airy disk model
init_amplitude = np.ptp(psf_subimage)
xcom, ycom = photutils.centroid_com(psf_subimage)
diam_1st_zero = (fwhm * 2.44) / 1.028
airy = models.AiryDisk2D(amplitude=init_amplitude, x_0=xcom, y_0=ycom,
radius=diam_1st_zero/2.)
# Levenberg-Marquardt algorithm
fitter = fitting.LevMarLSQFitter()
y, x = np.indices(psf_subimage.shape)
fit = fitter(airy, x, y, psf_subimage)
if crop:
mean_y = fit.y_0.value + suby
mean_x = fit.x_0.value + subx
else:
mean_y = fit.y_0.value
mean_x = fit.x_0.value
amplitude = fit.amplitude.value