diff --git a/src/autoprof/autoprofutils/SharedFunctions.py b/src/autoprof/autoprofutils/SharedFunctions.py index 9991c619..aa770801 100644 --- a/src/autoprof/autoprofutils/SharedFunctions.py +++ b/src/autoprof/autoprofutils/SharedFunctions.py @@ -473,6 +473,49 @@ def interpolate_bicubic(dat, X, Y): return f_interp(Y, X, grid=False) +def _cubic_lagrange_weights(t): + return np.array( + [ + -t * (t - 1) * (t - 2) / 6, + (t + 1) * (t - 1) * (t - 2) / 2, + -(t + 1) * t * (t - 2) / 2, + (t + 1) * t * (t - 1) / 6, + ] + ) + + +def _interpolate_bicubic_local(dat, X, Y, mask=None): + flux = np.full(len(X), np.nan) + valid = np.ones(len(X), dtype=bool) + + for i in range(len(X)): + if not np.isfinite(X[i]) or not np.isfinite(Y[i]): + valid[i] = False + continue + + xbase = int(np.floor(X[i])) + ybase = int(np.floor(Y[i])) + x0, x1 = xbase - 1, xbase + 3 + y0, y1 = ybase - 1, ybase + 3 + if x0 < 0 or y0 < 0 or x1 > dat.shape[1] or y1 > dat.shape[0]: + valid[i] = False + continue + + chunk = dat[y0:y1, x0:x1] + if not np.all(np.isfinite(chunk)): + valid[i] = False + continue + if mask is not None and np.any(mask[y0:y1, x0:x1]): + valid[i] = False + continue + + wx = _cubic_lagrange_weights(X[i] - xbase) + wy = _cubic_lagrange_weights(Y[i] - ybase) + flux[i] = np.sum(chunk * wy.reshape((-1, 1)) * wx) + + return flux, valid + + def interpolate_Lanczos(dat, X, Y, scale): """ Perform Lanczos interpolation on an image. @@ -517,6 +560,30 @@ def interpolate_Lanczos(dat, X, Y, scale): return np.array(flux) +def _interpolation_footprint_valid(dat, X, Y, mask=None, interp_method="lanczos", interp_window=5): + valid = np.ones(len(X), dtype=bool) + for i in range(len(X)): + if interp_method == "lanczos": + x0 = int(round(np.floor(X[i]) - interp_window + 1)) + x1 = int(round(np.floor(X[i]) + interp_window + 1)) + y0 = int(round(np.floor(Y[i]) - interp_window + 1)) + y1 = int(round(np.floor(Y[i]) + interp_window + 1)) + else: + raise ValueError( + "Unknown interpolate method %s. Should be lanczos" % interp_method + ) + if x0 < 0 or y0 < 0 or x1 > dat.shape[1] or y1 > dat.shape[0]: + valid[i] = False + continue + chunk = dat[y0:y1, x0:x1] + if chunk.size == 0 or not np.all(np.isfinite(chunk)): + valid[i] = False + continue + if mask is not None and np.any(mask[y0:y1, x0:x1]): + valid[i] = False + return valid + + def _iso_between( IMG, sma_low, @@ -560,38 +627,32 @@ def _iso_between( RR /= SuperEllipse_Rscale * Fmode_Rscale rselect = np.logical_and(RR < sma_high, RR > sma_low) fluxes = IMG[ranges[1][0] : ranges[1][1], ranges[0][0] : ranges[0][1]][rselect] - CHOOSE = None - if not mask is None and sma_high > 5: - CHOOSE = np.logical_not( - mask[ranges[1][0] : ranges[1][1], ranges[0][0] : ranges[0][1]][rselect] + theta_values = theta[rselect] + CHOOSE = np.isfinite(fluxes) + if not mask is None: + CHOOSE = np.logical_and( + CHOOSE, + np.logical_not( + mask[ranges[1][0] : ranges[1][1], ranges[0][0] : ranges[0][1]][rselect] + ), ) # Perform sigma clipping if requested if sigmaclip: - sclim = Sigma_Clip_Upper(fluxes, sclip_iterations, sclip_nsigma) - if CHOOSE is None: - CHOOSE = fluxes < sclim - else: - CHOOSE = np.logical_or(CHOOSE, fluxes < sclim) - if CHOOSE is not None and np.sum(CHOOSE) < 5: + if np.any(CHOOSE): + sclim = Sigma_Clip_Upper(fluxes[CHOOSE], sclip_iterations, sclip_nsigma) + clipped_choose = np.logical_and(CHOOSE, fluxes < sclim) + if np.any(clipped_choose): + CHOOSE = clipped_choose + countmasked = np.sum(np.logical_not(CHOOSE)) + if np.sum(CHOOSE) < 5: logging.warning( - "Entire Isophote is Masked! R_l: %.3f, R_h: %.3f, PA: %.3f, ellip: %.3f" + "Too few usable pixels in isophote band! R_l: %.3f, R_h: %.3f, PA: %.3f, ellip: %.3f" % (sma_low, sma_high, PARAMS["pa"] * 180 / np.pi, PARAMS["ellip"]) ) - CHOOSE = np.ones(CHOOSE.shape).astype(bool) - if CHOOSE is not None: - countmasked = np.sum(np.logical_not(CHOOSE)) - else: - countmasked = 0 if more: - if CHOOSE is not None and sma_high > 5: - return fluxes[CHOOSE], theta[rselect][CHOOSE], countmasked - else: - return fluxes, theta[rselect], countmasked + return fluxes[CHOOSE], theta_values[CHOOSE], countmasked else: - if CHOOSE is not None and sma_high > 5: - return fluxes[CHOOSE] - else: - return fluxes + return fluxes[CHOOSE] def _iso_extract( @@ -651,19 +712,20 @@ def _iso_extract( theta = theta[BORDER] Rlim = np.max(R) + footprint_choose = None if Rlim < rad_interp: - box = [ - [max(0, int(c["x"] - Rlim - 5)), min(IMG.shape[1], int(c["x"] + Rlim + 5))], - [max(0, int(c["y"] - Rlim - 5)), min(IMG.shape[0], int(c["y"] + Rlim + 5))], - ] if interp_method == "bicubic": - flux = interpolate_bicubic( - IMG[box[1][0] : box[1][1], box[0][0] : box[0][1]], - X - box[0][0], - Y - box[1][0], - ) + flux, footprint_choose = _interpolate_bicubic_local(IMG, X, Y, mask=mask) elif interp_method == "lanczos": flux = interpolate_Lanczos(IMG, X, Y, interp_window) + footprint_choose = _interpolation_footprint_valid( + IMG, + X, + Y, + mask=mask, + interp_method=interp_method, + interp_window=interp_window, + ) else: raise ValueError( "Unknown interpolate method %s. Should be one of lanczos or bicubic" % interp_method @@ -671,35 +733,43 @@ def _iso_extract( else: # round to integers and sample pixels values flux = IMG[np.rint(Y).astype(np.int32), np.rint(X).astype(np.int32)] - # CHOOSE holds bolean array for which flux values to keep, initialized as None for no clipping - CHOOSE = None + # CHOOSE holds bolean array for which flux values to keep. + # Interpolated samples are rejected when their interpolation footprint + # contains masked or non-finite pixels. + finite_flux = np.isfinite(flux) + CHOOSE = finite_flux.copy() # Mask pixels if a mask is given if not mask is None: - CHOOSE = np.logical_not(mask[np.rint(Y).astype(np.int32), np.rint(X).astype(np.int32)]) + CHOOSE = np.logical_and( + CHOOSE, np.logical_not(mask[np.rint(Y).astype(np.int32), np.rint(X).astype(np.int32)]) + ) + if footprint_choose is not None: + CHOOSE = np.logical_and(CHOOSE, footprint_choose) # Perform sigma clipping if requested if sigmaclip: - sclim = Sigma_Clip_Upper(flux, sclip_iterations, sclip_nsigma) - if CHOOSE is None: - CHOOSE = flux < sclim - else: - CHOOSE = np.logical_or(CHOOSE, flux < sclim) + if np.any(CHOOSE): + sclim = Sigma_Clip_Upper(flux[CHOOSE], sclip_iterations, sclip_nsigma) + clipped_choose = np.logical_and(CHOOSE, flux < sclim) + if np.any(clipped_choose): + CHOOSE = clipped_choose # Dont clip pixels if that removes all of the pixels - countmasked = 0 - if not CHOOSE is None and np.sum(CHOOSE) <= 0: + countmasked = np.sum(np.logical_not(CHOOSE)) + if np.sum(CHOOSE) <= 0: logging.warning( "Entire Isophote is Masked! R: %.3f, PA: %.3f, ellip: %.3f" % (sma, PARAMS["pa"] * 180 / np.pi, PARAMS["ellip"]) ) + flux = flux[CHOOSE] + theta = theta[CHOOSE] # Interpolate clipped flux values if requested - elif not CHOOSE is None and interp_mask: + elif interp_mask: flux[np.logical_not(CHOOSE)] = np.interp( theta[np.logical_not(CHOOSE)], theta[CHOOSE], flux[CHOOSE], period=2 * np.pi ) # simply remove all clipped pixels if user doesn't reqest another option - elif not CHOOSE is None: + elif not np.all(CHOOSE): flux = flux[CHOOSE] theta = theta[CHOOSE] - countmasked = np.sum(np.logical_not(CHOOSE)) # Return just the flux values, or flux and angle values if more: @@ -708,7 +778,7 @@ def _iso_extract( return flux -def _iso_line(IMG, length, width, pa, c, more=False): +def _iso_line(IMG, length, width, pa, c, more=False, mask=None): start = np.array([c["x"], c["y"]]) end = start + length * np.array([np.cos(pa), np.sin(pa)]) @@ -733,11 +803,21 @@ def _iso_line(IMG, length, width, pa, c, more=False): lselect = np.logical_and.reduce((XX >= -0.5, XX <= length, np.abs(YY) <= (width / 2))) flux = IMG[ranges[1][0] : ranges[1][1], ranges[0][0] : ranges[0][1]][lselect] + Xline = XX[lselect] + Yline = YY[lselect] + CHOOSE = np.isfinite(flux) + if mask is not None: + CHOOSE = np.logical_and( + CHOOSE, + np.logical_not( + mask[ranges[1][0] : ranges[1][1], ranges[0][0] : ranges[0][1]][lselect] + ), + ) if more: - return flux, XX[lselect], YY[lselect] + return flux[CHOOSE], Xline[CHOOSE], Yline[CHOOSE] else: - return flux, XX[lselect] + return flux[CHOOSE], Xline[CHOOSE] def StarFind( @@ -970,6 +1050,15 @@ def StarFind( } +def _photutils_masked_data(dat, results): + mask = np.logical_not(np.isfinite(dat)) + if results.get("mask", None) is not None: + mask = np.logical_or(mask, results["mask"]) + if np.any(mask): + return np.ma.array(dat, mask=mask) + return dat + + def _x_to_pa(x): """ Internal, basic function to ensure position angles remain diff --git a/src/autoprof/pipeline_steps/Axial_Profiles.py b/src/autoprof/pipeline_steps/Axial_Profiles.py index df1ee434..bd54489e 100644 --- a/src/autoprof/pipeline_steps/Axial_Profiles.py +++ b/src/autoprof/pipeline_steps/Axial_Profiles.py @@ -145,6 +145,7 @@ def Axial_Profiles(IMG, results, options): "y": results["center"]["y"] + ang * rd * pR * np.sin(pa + (0 if ang > 0 else np.pi)), }, + mask=mask, ) for oi, oR in enumerate(R): length = (R[oi] - R[oi - 1]) if oi > 0 else 1.0 diff --git a/src/autoprof/pipeline_steps/Background.py b/src/autoprof/pipeline_steps/Background.py index f79a94ab..c606baa2 100644 --- a/src/autoprof/pipeline_steps/Background.py +++ b/src/autoprof/pipeline_steps/Background.py @@ -1,8 +1,12 @@ -from photutils.segmentation import SegmentationImage +from photutils.segmentation import make_source_mask +from astropy.utils.exceptions import AstropyDeprecationWarning from scipy.stats import iqr from scipy.fftpack import fft2, ifft2 +from scipy.interpolate import SmoothBivariateSpline +from scipy.ndimage import distance_transform_edt, label import logging import numpy as np +import warnings from ..autoprofutils.SharedFunctions import Smooth_Mode from ..autoprofutils.Diagnostic_Plots import Plot_Background @@ -16,6 +20,101 @@ ) +def _background_border_mask(shape): + mask = np.ones(shape, dtype=bool) + mask[ + int(shape[0] / 5.0) : int(4.0 * shape[0] / 5.0), + int(shape[1] / 5.0) : int(4.0 * shape[1] / 5.0), + ] = False + return mask + + +def _finite_sample_values(IMG, sample_mask): + values = IMG[sample_mask].flatten() + return values[np.isfinite(values)] + + +def _require_background_values(values): + if len(values) == 0: + raise ValueError("No finite pixels available for background estimation") + return values + + +def _apply_background_speedup(values, options): + if "ap_background_speedup" in options and int(options["ap_background_speedup"]) > 1: + values = values[:: int(options["ap_background_speedup"])] + return values + + +def _nearest_finite_fill(IMG, bad): + filled = np.array(IMG, dtype=float, copy=True) + if np.all(bad): + raise ValueError("No finite pixels available for background estimation") + dumy, indices = distance_transform_edt( + bad, + return_distances=True, + return_indices=True, + ) + filled[bad] = filled[tuple(indices[:, bad])] + return filled + + +def _spline_fill_nonfinite(IMG, options): + bad = np.logical_not(np.isfinite(IMG)) + if not np.any(bad): + return IMG + + labels, nlabels = label(bad) + if nlabels > 0: + largest_bad = np.max(np.bincount(labels.ravel())[1:]) + else: + largest_bad = 0 + bad_fraction = np.sum(bad) / bad.size + if bad_fraction > 0.01 or largest_bad > 0.005 * bad.size: + logging.warning( + "%s: unsharp background has large non-finite regions; spline-filled FFT background may be unreliable" + % options["ap_name"] + ) + + filled = np.array(IMG, dtype=float, copy=True) + good = np.logical_not(bad) + y, x = np.nonzero(good) + z = filled[good] + y_bad, x_bad = np.nonzero(bad) + + kx = min(3, len(np.unique(x)) - 1) + ky = min(3, len(np.unique(y)) - 1) + if kx < 1 or ky < 1: + return _nearest_finite_fill(IMG, bad) + + max_spline_points = 50000 + if len(z) > max_spline_points: + step = int(np.ceil(len(z) / max_spline_points)) + x_fit = x[::step] + y_fit = y[::step] + z_fit = z[::step] + else: + x_fit = x + y_fit = y + z_fit = z + + try: + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + spline = SmoothBivariateSpline(x_fit, y_fit, z_fit, kx=kx, ky=ky) + fill_values = spline.ev(x_bad, y_bad) + if not np.all(np.isfinite(fill_values)): + raise ValueError("Spline produced non-finite fill values") + filled[bad] = fill_values + return filled + except Exception as e: + logging.warning( + "%s: spline fill failed for unsharp background (%s); using nearest finite fill" + % (options["ap_name"], str(e)) + ) + return _nearest_finite_fill(IMG, bad) + + def Background_Mode(IMG, results, options): """Compute the mode flux in the border of an image. @@ -64,24 +163,19 @@ def Background_Mode(IMG, results, options): """ # Mask main body of image so only outer 1/5th is used # for background calculation. - if "mask" in results and not results["mask"] is None and np.any(results["mask"]): + use_existing_mask = "mask" in results and not results["mask"] is None and np.any(results["mask"]) + if use_existing_mask: mask = np.logical_not(results["mask"]) logging.info( "%s: Background using mask. Masking %i pixels" % (options["ap_name"], np.sum(results["mask"])) ) else: - mask = np.ones(IMG.shape, dtype=bool) - mask[ - int(IMG.shape[0] / 5.0) : int(4.0 * IMG.shape[0] / 5.0), - int(IMG.shape[1] / 5.0) : int(4.0 * IMG.shape[1] / 5.0), - ] = False - values = IMG[mask].flatten() - if len(values) < 1e5: - values = IMG.flatten() - if "ap_background_speedup" in options and int(options["ap_background_speedup"]) > 1: - values = values[:: int(options["ap_background_speedup"])] - values = values[np.isfinite(values)] + mask = _background_border_mask(IMG.shape) + values = _finite_sample_values(IMG, mask) + if len(values) < 1e5 and not use_existing_mask: + values = IMG[np.isfinite(IMG)].flatten() + values = _require_background_values(_apply_background_speedup(values, options)) if "ap_set_background" in options: bkgrnd = options["ap_set_background"] @@ -95,7 +189,12 @@ def Background_Mode(IMG, results, options): noise = options["ap_set_background_noise"] logging.info("%s: Background Noise set by user: %.4e" % (options["ap_name"], noise)) else: - noise = iqr(values[(values - bkgrnd) < 0], rng=[100 - 68.2689492137, 100]) + noise_values = values[(values - bkgrnd) < 0] + noise = ( + iqr(noise_values, rng=[100 - 68.2689492137, 100]) + if len(noise_values) > 0 + else np.nan + ) if not np.isfinite(noise): noise = iqr(values, rng=[16, 84]) / 2.0 uncertainty = noise / np.sqrt(np.sum((values - bkgrnd) < 0)) @@ -150,44 +249,55 @@ def Background_DilatedSources(IMG, results, options): # Mask main body of image so only outer 1/5th is used # for background calculation. if "mask" in results and not results["mask"] is None: - mask = results["mask"] + mask = results["mask"].astype(bool) else: - mask = np.zeros(IMG.shape) + mask = np.zeros(IMG.shape, dtype=bool) mask[ int(IMG.shape[0] / 5.0) : int(4.0 * IMG.shape[0] / 5.0), int(IMG.shape[1] / 5.0) : int(4.0 * IMG.shape[1] / 5.0), - ] = 1 + ] = True + mask = np.logical_or(mask, np.logical_not(np.isfinite(IMG))) + base_values = _require_background_values(_finite_sample_values(IMG, np.logical_not(mask))) # Run photutils source mask to remove pixels with sources # such as stars and galaxies, including a boarder # around each source. if not ("ap_set_background" in options and "ap_set_background_noise" in options): - segm = SegmentationImage(IMG) - source_mask = segm.make_source_mask( - nsigma=3, - npixels=int(1.0 / options["ap_pixscale"]), - dilate_size=40, - filter_fwhm=1.0 / options["ap_pixscale"], - filter_size=int(3.0 / options["ap_pixscale"]), - sigclip_iters=5, - ) + finite_values = _require_background_values(IMG[np.isfinite(IMG)].flatten()) + finite_img = np.array(IMG, copy=True) + finite_img[np.logical_not(np.isfinite(finite_img))] = np.median(finite_values) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", AstropyDeprecationWarning) + source_mask = make_source_mask( + finite_img, + nsigma=3, + npixels=int(1.0 / options["ap_pixscale"]), + mask=mask, + dilate_size=40, + filter_fwhm=1.0 / options["ap_pixscale"], + filter_size=int(3.0 / options["ap_pixscale"]), + sigclip_iters=5, + ) mask = np.logical_or(mask, source_mask) + values = _finite_sample_values(IMG, np.logical_not(mask)) + if len(values) == 0: + values = base_values # Return statistics from background sky bkgrnd = ( options["ap_set_background"] if "ap_set_background" in options - else np.median(IMG[np.logical_not(mask)]) + else np.median(values) ) noise = ( options["ap_set_background_noise"] if "ap_set_background_noise" in options - else iqr(IMG[np.logical_not(mask)], rng=[16, 84]) / 2 + else iqr(values, rng=[16, 84]) / 2 ) - uncertainty = noise / np.sqrt(np.sum(np.logical_not(mask))) + uncertainty = noise / np.sqrt(len(values)) if "ap_doplot" in options and options["ap_doplot"]: - Plot_Background(IMG[np.logical_not(mask)].ravel(), bkgrnd, noise, results, options) + Plot_Background(values, bkgrnd, noise, results, options) return IMG, { "background": bkgrnd, "background noise": noise, @@ -241,12 +351,10 @@ def Background_Basic(IMG, results, options): int(IMG.shape[0] / 5.0) : int(4.0 * IMG.shape[0] / 5.0), int(IMG.shape[1] / 5.0) : int(4.0 * IMG.shape[1] / 5.0), ] = False - values = IMG[mask].flatten() + values = _finite_sample_values(IMG, mask) if len(values) < 1e3: - values = IMG.flatten() - if "ap_background_speedup" in options and int(options["ap_background_speedup"]) > 1: - values = values[:: int(options["ap_background_speedup"])] - values = values[np.isfinite(values)] + values = IMG[np.isfinite(IMG)].flatten() + values = _require_background_values(_apply_background_speedup(values, options)) bkgrnd = options["ap_set_background"] if "ap_set_background" in options else np.mean(values) noise = ( @@ -297,7 +405,9 @@ def Background_Unsharp(IMG, results, options): """ - coefs = fft2(IMG) + dumy, stats = Background_Mode(IMG, results, options) + finite_img = _spline_fill_nonfinite(IMG, options) + coefs = fft2(finite_img) unsharp = ( int(options["ap_background_unsharp_lowpass"]) @@ -307,6 +417,5 @@ def Background_Unsharp(IMG, results, options): coefs[unsharp:-unsharp] = 0 coefs[:, unsharp:-unsharp] = 0 - dumy, stats = Background_Mode(IMG, results, options) stats.update({"background": ifft2(coefs).real}) return IMG, stats diff --git a/src/autoprof/pipeline_steps/Center.py b/src/autoprof/pipeline_steps/Center.py index b6923d1a..7fb03c41 100644 --- a/src/autoprof/pipeline_steps/Center.py +++ b/src/autoprof/pipeline_steps/Center.py @@ -20,6 +20,16 @@ __all__ = ("Center_Forced", "Center_2DGaussian", "Center_1DGaussian", "Center_OfMass", "Center_Peak", "Center_HillClimb", "Center_HillClimb_mean") + +def _center_mask(dat, results, extra_mask=None): + center_mask = np.logical_not(np.isfinite(dat)) + if results.get("mask", None) is not None: + center_mask = np.logical_or(center_mask, results["mask"]) + if extra_mask is not None: + center_mask = np.logical_or(center_mask, extra_mask) + return center_mask + + def Center_Forced(IMG, results, options): """Extracts previously fit center coordinates. @@ -81,6 +91,7 @@ def Center_Forced(IMG, results, options): """ current_center = {"x": IMG.shape[1] / 2, "y": IMG.shape[0] / 2} + dat = IMG - results["background"] if "ap_guess_center" in options: current_center = deepcopy(options["ap_guess_center"]) logging.info( @@ -92,16 +103,7 @@ def Center_Forced(IMG, results, options): "%s: Center set by user: %s" % (options["ap_name"], str(options["ap_set_center"])) ) - sb0 = flux_to_sb( - _iso_extract( - IMG - results["background"], - 0.0, - {"ellip": 0.0, "pa": 0.0}, - options["ap_set_center"], - )[0], - options["ap_pixscale"], - options["ap_zeropoint"] if "ap_zeropoint" in options else 22.5, - ) + sb0 = _central_surface_brightness(dat, options["ap_set_center"], results, options) return IMG, { "center": deepcopy(options["ap_set_center"]), "auxfile central sb": "central surface brightness: %.4f mag arcsec^-2" @@ -129,13 +131,7 @@ def Center_Forced(IMG, results, options): "%s: Forced center failed! Using image center (or guess)." % options["ap_name"] ) - sb0 = flux_to_sb( - _iso_extract( - IMG - results["background"], 0.0, {"ellip": 0.0, "pa": 0.0}, current_center - )[0], - options["ap_pixscale"], - options["ap_zeropoint"] if "ap_zeropoint" in options else 22.5, - ) + sb0 = _central_surface_brightness(dat, current_center, results, options) return IMG, { "center": current_center, "auxfile center": "center x: %.2f pix, y: %.2f pix" @@ -219,6 +215,8 @@ def Center_2DGaussian(IMG, results, options): ) return IMG, {"center": deepcopy(options["ap_set_center"])} + dat = IMG - results["background"] + # Create mask to focus centering algorithm on the center of the image ranges = [ [ @@ -276,11 +274,16 @@ def Center_2DGaussian(IMG, results, options): ] centralize_mask = np.ones(IMG.shape, dtype=bool) centralize_mask[ranges[1][0] : ranges[1][1], ranges[0][0] : ranges[0][1]] = False + center_mask = _center_mask(dat, results, centralize_mask) try: - x, y = centroid_2dg(IMG - results["background"], mask=centralize_mask) + if not np.any(np.logical_not(center_mask)): + raise ValueError("all center pixels are masked") + x, y = centroid_2dg(np.where(center_mask, 0.0, dat), mask=center_mask) + if not np.all(np.isfinite([x, y])): + raise ValueError("center is non-finite") current_center = {"x": x, "y": y} - except: + except Exception: logging.warning( "%s: 2D Gaussian center finding failed! using image center (or guess)." % options["ap_name"] @@ -289,7 +292,7 @@ def Center_2DGaussian(IMG, results, options): # Plot center value for diagnostic purposes if "ap_doplot" in options and options["ap_doplot"]: plt.imshow( - np.clip(IMG - results["background"], a_min=0, a_max=None), + np.clip(dat, a_min=0, a_max=None), origin="lower", cmap="Greys_r", norm=ImageNormalize(stretch=LogStretch()), @@ -306,7 +309,10 @@ def Center_2DGaussian(IMG, results, options): dpi=options["ap_plotdpi"] if "ap_plotdpi" in options else 300, ) plt.close() - logging.info("%s Center found: x %.1f, y %.1f" % (options["ap_name"], x, y)) + logging.info( + "%s Center found: x %.1f, y %.1f" + % (options["ap_name"], current_center["x"], current_center["y"]) + ) return IMG, { "center": current_center, "auxfile center": "center x: %.2f pix, y: %.2f pix" @@ -391,6 +397,8 @@ def Center_1DGaussian(IMG, results, options): ) return IMG, {"center": deepcopy(options["ap_set_center"])} + dat = IMG - results["background"] + # Create mask to focus centering algorithm on the center of the image ranges = [ [ @@ -448,31 +456,45 @@ def Center_1DGaussian(IMG, results, options): ] centralize_mask = np.ones(IMG.shape, dtype=bool) centralize_mask[ranges[1][0] : ranges[1][1], ranges[0][0] : ranges[0][1]] = False + center_mask = _center_mask(dat, results, centralize_mask) try: - x, y = centroid_1dg(IMG - results["background"], mask=centralize_mask) + if not np.any(np.logical_not(center_mask)): + raise ValueError("all center pixels are masked") + x, y = centroid_1dg(np.where(center_mask, 0.0, dat), mask=center_mask) + if not np.all(np.isfinite([x, y])): + raise ValueError("center is non-finite") current_center = {"x": x, "y": y} - except: + except Exception: logging.warning( - "%s: 2D Gaussian center finding failed! using image center (or guess)." + "%s: 1D Gaussian center finding failed! using image center (or guess)." % options["ap_name"] ) # Plot center value for diagnostic purposes if "ap_doplot" in options and options["ap_doplot"]: plt.imshow( - np.clip(IMG - results["background"], a_min=0, a_max=None), + np.clip(dat, a_min=0, a_max=None), origin="lower", cmap="Greys_r", norm=ImageNormalize(stretch=LogStretch()), ) - plt.plot([y], [x], marker="x", markersize=10, color="y") + plt.plot( + [current_center["x"]], + [current_center["y"]], + marker="x", + markersize=10, + color="y", + ) plt.savefig( f"{options.get('ap_plotpath','')}center_vis_{options['ap_name']}.{options.get('ap_plot_extension', 'jpg')}", dpi=options["ap_plotdpi"] if "ap_plotdpi" in options else 300, ) plt.close() - logging.info("%s Center found: x %.1f, y %.1f" % (options["ap_name"], x, y)) + logging.info( + "%s Center found: x %.1f, y %.1f" + % (options["ap_name"], current_center["x"], current_center["y"]) + ) return IMG, { "center": current_center, "auxfile center": "center x: %.2f pix, y: %.2f pix" @@ -558,13 +580,7 @@ def Center_OfMass(IMG, results, options): "%s: Center set by user: %s" % (options["ap_name"], str(options["ap_set_center"])) ) - sb0 = flux_to_sb( - _iso_extract(dat, 0.0, {"ellip": 0.0, "pa": 0.0}, options["ap_set_center"])[ - 0 - ], - options["ap_pixscale"], - options["ap_zeropoint"] if "ap_zeropoint" in options else 22.5, - ) + sb0 = _central_surface_brightness(dat, options["ap_set_center"], results, options) return IMG, { "center": deepcopy(options["ap_set_center"]), "auxfile central sb": "central surface brightness: %.4f mag arcsec^-2" @@ -575,7 +591,7 @@ def Center_OfMass(IMG, results, options): (options["ap_centeringring"] if "ap_centeringring" in options else 10) * results["psf fwhm"] ) - xx, yy = np.meshgrid(np.arange(searchring), np.arange(searchring)) + center_mask = _center_mask(dat, results) N_updates = 0 while N_updates < 100: N_updates += 1 @@ -590,29 +606,46 @@ def Center_OfMass(IMG, results, options): min(IMG.shape[0], int(current_center["y"] + searchring / 2)), ], ] - current_center = { - "x": ranges[0][0] - + np.sum( - (dat[ranges[1][0] : ranges[1][1], ranges[0][0] : ranges[0][1]] * xx) + chunk = dat[ranges[1][0] : ranges[1][1], ranges[0][0] : ranges[0][1]] + chunk_mask = center_mask[ + ranges[1][0] : ranges[1][1], ranges[0][0] : ranges[0][1] + ] + choose = np.logical_not(chunk_mask) + if not np.any(choose): + logging.warning( + "%s: Center of mass failed because all center pixels are masked." + % options["ap_name"] ) - / np.sum(dat[ranges[1][0] : ranges[1][1], ranges[0][0] : ranges[0][1]]), - "y": ranges[1][0] - + np.sum( - (dat[ranges[1][0] : ranges[1][1], ranges[0][0] : ranges[0][1]] * yy) + break + yy, xx = np.indices(chunk.shape) + weights = np.where(choose, chunk, 0.0) + denominator = np.sum(weights) + if (not np.isfinite(denominator)) or denominator == 0: + logging.warning( + "%s: Center of mass failed because the center flux sum is invalid." + % options["ap_name"] ) - / np.sum(dat[ranges[1][0] : ranges[1][1], ranges[0][0] : ranges[0][1]]), + current_center = old_center + break + new_center = { + "x": ranges[0][0] + np.sum(weights * xx) / denominator, + "y": ranges[1][0] + np.sum(weights * yy) / denominator, } + if not np.all(np.isfinite([new_center["x"], new_center["y"]])): + logging.warning( + "%s: Center of mass failed because the center is non-finite." + % options["ap_name"] + ) + current_center = old_center + break + current_center = new_center if ( abs(current_center["x"] - old_center["x"]) < 0.1 * results["psf fwhm"] and abs(current_center["y"] - old_center["y"]) < 0.1 * results["psf fwhm"] ): break - sb0 = flux_to_sb( - _iso_extract(dat, 0.0, {"ellip": 0.0, "pa": 0.0}, current_center)[0], - options["ap_pixscale"], - options["ap_zeropoint"] if "ap_zeropoint" in options else 22.5, - ) + sb0 = _central_surface_brightness(dat, current_center, results, options) return IMG, { "center": current_center, "auxfile center": "center x: %.2f pix, y: %.2f pix" @@ -636,13 +669,7 @@ def Center_Peak(IMG, results, options): "%s: Center set by user: %s" % (options["ap_name"], str(options["ap_set_center"])) ) - sb0 = flux_to_sb( - _iso_extract(dat, 0.0, {"ellip": 0.0, "pa": 0.0}, options["ap_set_center"])[ - 0 - ], - options["ap_pixscale"], - options["ap_zeropoint"] if "ap_zeropoint" in options else 22.5, - ) + sb0 = _central_surface_brightness(dat, options["ap_set_center"], results, options) return IMG, { "center": deepcopy(options["ap_set_center"]), "auxfile central sb": "central surface brightness: %.4f mag arcsec^-2" @@ -653,22 +680,6 @@ def Center_Peak(IMG, results, options): (options["ap_centeringring"] if "ap_centeringring" in options else 10) * results["psf fwhm"] ) - xx, yy = np.meshgrid(np.arange(searchring), np.arange(searchring)) - xx = xx.flatten() - yy = yy.flatten() - A = np.array( - [ - np.ones(xx.shape), - xx, - yy, - xx ** 2, - yy ** 2, - xx * yy, - xx * yy ** 2, - yy * xx ** 2, - xx ** 2 * yy ** 2, - ] - ).T ranges = [ [ max(0, int(current_center["x"] - searchring / 2)), @@ -679,23 +690,49 @@ def Center_Peak(IMG, results, options): min(IMG.shape[0], int(current_center["y"] + searchring / 2)), ], ] - chunk = np.clip( - dat[ranges[1][0] : ranges[1][1], ranges[0][0] : ranges[0][1]].T, - a_min=results["background noise"] / 5, - a_max=None, - ) + center_mask = _center_mask(dat, results) + chunk = dat[ranges[1][0] : ranges[1][1], ranges[0][0] : ranges[0][1]] + chunk_mask = center_mask[ + ranges[1][0] : ranges[1][1], ranges[0][0] : ranges[0][1] + ] + yy, xx = np.indices(chunk.shape) + xx = xx.flatten() + yy = yy.flatten() + choose = np.logical_not(chunk_mask.flatten()) + floor = max(results["background noise"] / 5, np.finfo(float).tiny) + flux = np.clip(chunk.flatten()[choose], a_min=floor, a_max=None) - poly2dfit = np.linalg.lstsq(A, np.log10(chunk.flatten()), rcond=None) - current_center = { - "x": -poly2dfit[0][2] / (2 * poly2dfit[0][4]) + ranges[0][0], - "y": -poly2dfit[0][1] / (2 * poly2dfit[0][3]) + ranges[1][0], - } + try: + if len(flux) < 9: + raise ValueError("not enough unmasked pixels") + A = np.array( + [ + np.ones(xx[choose].shape), + xx[choose], + yy[choose], + xx[choose] ** 2, + yy[choose] ** 2, + xx[choose] * yy[choose], + xx[choose] * yy[choose] ** 2, + yy[choose] * xx[choose] ** 2, + xx[choose] ** 2 * yy[choose] ** 2, + ] + ).T + poly2dfit = np.linalg.lstsq(A, np.log10(flux), rcond=None) + new_center = { + "x": -poly2dfit[0][1] / (2 * poly2dfit[0][3]) + ranges[0][0], + "y": -poly2dfit[0][2] / (2 * poly2dfit[0][4]) + ranges[1][0], + } + if not np.all(np.isfinite([new_center["x"], new_center["y"]])): + raise ValueError("center is non-finite") + current_center = new_center + except Exception: + logging.warning( + "%s: Peak center finding failed! using image center (or guess)." + % options["ap_name"] + ) - sb0 = flux_to_sb( - _iso_extract(dat, 0.0, {"ellip": 0.0, "pa": 0.0}, current_center)[0], - options["ap_pixscale"], - options["ap_zeropoint"] if "ap_zeropoint" in options else 22.5, - ) + sb0 = _central_surface_brightness(dat, current_center, results, options) return IMG, { "center": current_center, "auxfile center": "center x: %.2f pix, y: %.2f pix" @@ -704,7 +741,24 @@ def Center_Peak(IMG, results, options): } -def _hillclimb_loss(x, IMG, PSF, noise): +def _central_surface_brightness(dat, center, results, options): + isovals = _iso_extract( + dat, + 0.0, + {"ellip": 0.0, "pa": 0.0}, + center, + mask=results.get("mask", None), + ) + if len(isovals) == 0: + return np.nan + return flux_to_sb( + isovals[0], + options["ap_pixscale"], + options["ap_zeropoint"] if "ap_zeropoint" in options else 22.5, + ) + + +def _hillclimb_loss(x, IMG, PSF, noise, mask=None): center_loss = 0 for rr in range(3): RR = (rr + 1.0) * PSF / 2 @@ -724,11 +778,17 @@ def _hillclimb_loss(x, IMG, PSF, noise): rad_interp=10 * PSF, interp_method="lanczos", interp_window=3, + mask=mask, ) + if len(isovals) < 2: + return np.inf coefs = fft(isovals) - center_loss += np.abs(coefs[1]) / ( - len(isovals) * (max(0, np.median(isovals)) + noise) - ) + denominator = len(isovals) * (max(0, np.median(isovals)) + noise) + if (not np.isfinite(denominator)) or denominator <= 0: + return np.inf + center_loss += np.abs(coefs[1]) / denominator + if not np.isfinite(center_loss): + return np.inf return center_loss @@ -807,13 +867,7 @@ def Center_HillClimb(IMG, results, options): "%s: Center set by user: %s" % (options["ap_name"], str(options["ap_set_center"])) ) - sb0 = flux_to_sb( - _iso_extract(dat, 0.0, {"ellip": 0.0, "pa": 0.0}, options["ap_set_center"], mask = results.get("mask", None))[ - 0 - ], - options["ap_pixscale"], - options["ap_zeropoint"] if "ap_zeropoint" in options else 22.5, - ) + sb0 = _central_surface_brightness(dat, options["ap_set_center"], results, options) return IMG, { "center": deepcopy(options["ap_set_center"]), "auxfile central sb": "central surface brightness: %.4f mag arcsec^-2" @@ -828,18 +882,25 @@ def Center_HillClimb(IMG, results, options): track_centers = [] small_update_count = 0 total_count = 0 + refine_center = True while small_update_count <= 5 and total_count <= 100: total_count += 1 phases = [] isovals = [] coefs = [] + sampled_radii = [] for r in sampleradii: - isovals.append( - _iso_extract( - dat, r, {"ellip": 0.0, "pa": 0.0}, current_center, more=True, - mask = results.get("mask", None) - ) + isovals_r = _iso_extract( + dat, + r, + {"ellip": 0.0, "pa": 0.0}, + current_center, + more=True, + mask=results.get("mask", None), ) + if len(isovals_r[0]) < 2: + continue + isovals.append(isovals_r) coefs.append( fft( np.clip( @@ -850,10 +911,25 @@ def Center_HillClimb(IMG, results, options): ) ) phases.append((-np.angle(coefs[-1][1])) % (2 * np.pi)) + sampled_radii.append(r) + if len(phases) == 0: + logging.warning( + "%s: Center finding stopped because all sampled isophotes were masked" + % options["ap_name"] + ) + refine_center = False + break direction = Angle_Median(phases) % (2 * np.pi) + if not np.isfinite(direction): + logging.warning( + "%s: Center finding stopped because the update direction is invalid" + % options["ap_name"] + ) + refine_center = False + break levels = [] level_locs = [] - for i, r in enumerate(sampleradii): + for i, r in enumerate(sampled_radii): floc = np.argmin(np.abs((isovals[i][1] % (2 * np.pi)) - direction)) rloc = np.argmin( np.abs( @@ -883,38 +959,41 @@ def Center_HillClimb(IMG, results, options): small_update_count = 0 track_centers.append([current_center["x"], current_center["y"]]) - # refine center - ranges = [ - [ - max(0, int(current_center["x"] - results["psf fwhm"] * 5)), - min(dat.shape[1], int(current_center["x"] + results["psf fwhm"] * 5)), - ], - [ - max(0, int(current_center["y"] - results["psf fwhm"] * 5)), - min(dat.shape[0], int(current_center["y"] + results["psf fwhm"] * 5)), - ], - ] - - res = minimize( - _hillclimb_loss, - x0=[current_center["x"] - ranges[0][0], current_center["y"] - ranges[1][0]], - args=( - dat[ranges[1][0] : ranges[1][1], ranges[0][0] : ranges[0][1]], - results["psf fwhm"], - results["background noise"], - ), - method="Nelder-Mead", - ) - if res.success: - current_center["x"] = res.x[0] + ranges[0][0] - current_center["y"] = res.x[1] + ranges[1][0] + if refine_center: + # refine center + ranges = [ + [ + max(0, int(current_center["x"] - results["psf fwhm"] * 5)), + min(dat.shape[1], int(current_center["x"] + results["psf fwhm"] * 5)), + ], + [ + max(0, int(current_center["y"] - results["psf fwhm"] * 5)), + min(dat.shape[0], int(current_center["y"] + results["psf fwhm"] * 5)), + ], + ] + refine_mask = None + if "mask" in results: + refine_mask = results["mask"][ + ranges[1][0] : ranges[1][1], ranges[0][0] : ranges[0][1] + ] + + res = minimize( + _hillclimb_loss, + x0=[current_center["x"] - ranges[0][0], current_center["y"] - ranges[1][0]], + args=( + dat[ranges[1][0] : ranges[1][1], ranges[0][0] : ranges[0][1]], + results["psf fwhm"], + results["background noise"], + refine_mask, + ), + method="Nelder-Mead", + ) + if res.success and np.all(np.isfinite(res.x)) and np.isfinite(res.fun): + current_center["x"] = res.x[0] + ranges[0][0] + current_center["y"] = res.x[1] + ranges[1][0] track_centers.append([current_center["x"], current_center["y"]]) - sb0 = flux_to_sb( - _iso_extract(dat, 0.0, {"ellip": 0.0, "pa": 0.0}, current_center, mask = results.get("mask", None))[0], - options["ap_pixscale"], - options["ap_zeropoint"] if "ap_zeropoint" in options else 22.5, - ) + sb0 = _central_surface_brightness(dat, current_center, results, options) return IMG, { "center": current_center, "auxfile center": "center x: %.2f pix, y: %.2f pix" @@ -923,7 +1002,7 @@ def Center_HillClimb(IMG, results, options): } -def _hillclimb_mean_loss(x, IMG, PSF, noise): +def _hillclimb_mean_loss(x, IMG, PSF, noise, mask=None): center_loss = 0 for rr in range(3): isovals = _iso_extract( @@ -933,11 +1012,17 @@ def _hillclimb_mean_loss(x, IMG, PSF, noise): {"x": x[0], "y": x[1]}, more=False, rad_interp=10 * PSF, + mask=mask, ) + if len(isovals) < 2: + return np.inf coefs = fft(isovals) - center_loss += np.abs(coefs[1]) / ( - len(isovals) * (max(noise, np.mean(isovals))) - ) + denominator = len(isovals) * max(noise, np.mean(isovals)) + if (not np.isfinite(denominator)) or denominator <= 0: + return np.inf + center_loss += np.abs(coefs[1]) / denominator + if not np.isfinite(center_loss): + return np.inf return center_loss @@ -1031,23 +1116,46 @@ def Center_HillClimb_mean(IMG, results, options): track_centers = [] small_update_count = 0 total_count = 0 + refine_center = True while small_update_count <= 5 and total_count <= 100: total_count += 1 phases = [] isovals = [] coefs = [] + sampled_radii = [] for r in sampleradii: - isovals.append( - _iso_extract( - dat, r, {"ellip": 0.0, "pa": 0.0}, current_center, more=True - ) + isovals_r = _iso_extract( + dat, + r, + {"ellip": 0.0, "pa": 0.0}, + current_center, + more=True, + mask=results.get("mask", None), ) + if len(isovals_r[0]) < 2: + continue + isovals.append(isovals_r) coefs.append(fft(isovals[-1][0])) phases.append((-np.angle(coefs[-1][1])) % (2 * np.pi)) + sampled_radii.append(r) + if len(phases) == 0: + logging.warning( + "%s: Mean center finding stopped because all sampled isophotes were masked" + % options["ap_name"] + ) + refine_center = False + break direction = Angle_Median(phases) % (2 * np.pi) + if not np.isfinite(direction): + logging.warning( + "%s: Mean center finding stopped because the update direction is invalid" + % options["ap_name"] + ) + refine_center = False + break levels = [] level_locs = [] - for i, r in enumerate(sampleradii): + for i, r in enumerate(sampled_radii): floc = np.argmin(np.abs(isovals[i][1] - direction)) rloc = np.argmin( np.abs(isovals[i][1] - ((direction + np.pi) % (2 * np.pi))) @@ -1075,16 +1183,22 @@ def Center_HillClimb_mean(IMG, results, options): small_update_count = 0 track_centers.append([current_center["x"], current_center["y"]]) - # refine center - res = minimize( - _hillclimb_mean_loss, - x0=[current_center["x"], current_center["y"]], - args=(dat, results["psf fwhm"], results["background noise"]), - method="Nelder-Mead", - ) - if res.success: - current_center["x"] = res.x[0] - current_center["y"] = res.x[1] + if refine_center: + # refine center + res = minimize( + _hillclimb_mean_loss, + x0=[current_center["x"], current_center["y"]], + args=( + dat, + results["psf fwhm"], + results["background noise"], + results.get("mask", None), + ), + method="Nelder-Mead", + ) + if res.success and np.all(np.isfinite(res.x)) and np.isfinite(res.fun): + current_center["x"] = res.x[0] + current_center["y"] = res.x[1] return IMG, { "center": current_center, diff --git a/src/autoprof/pipeline_steps/Check_Fit.py b/src/autoprof/pipeline_steps/Check_Fit.py index 2fa42345..21368403 100644 --- a/src/autoprof/pipeline_steps/Check_Fit.py +++ b/src/autoprof/pipeline_steps/Check_Fit.py @@ -83,12 +83,18 @@ def Check_Fit(IMG, results, options): tests = {} # subtract background from image during processing dat = IMG - results["background"] + mask = results["mask"] if "mask" in results else None + if mask is not None and not np.any(mask): + mask = None # Compare variability of flux values along isophotes ###################################################################### use_center = results["center"] count_variable = 0 count_initrelative = 0 + count_checked = 0 + count_initrelative_checked = 0 + count_skipped = 0 f2_compare = [] f1_compare = [] if "fit R" in results: @@ -109,46 +115,79 @@ def Check_Fit(IMG, results, options): dat, checkson["R"][i], { - "ellip": results["init ellip"], # fixme, use mask + "ellip": results["init ellip"], "pa": results["init pa"], }, use_center, + mask=mask, ) + init_isovals = init_isovals[np.isfinite(init_isovals)] isovals = _iso_extract( dat, checkson["R"][i], {"ellip": checkson["ellip"][i], "pa": checkson["pa"][i]}, use_center, + mask=mask, ) + isovals = isovals[np.isfinite(isovals)] + if len(isovals) <= 2: + count_skipped += 1 + continue + med_isovals = np.median(isovals) + iqr_isovals = iqr(isovals) + if not (np.isfinite(med_isovals) and np.isfinite(iqr_isovals)): + count_skipped += 1 + continue coefs = fft(np.clip(isovals, a_max=np.quantile(isovals, 0.85), a_min=None)) + count_checked += 1 - if np.median(isovals) < (iqr(isovals) - results["background noise"]): + if med_isovals < (iqr_isovals - results["background noise"]): count_variable += 1 - if ( - (iqr(isovals) - results["background noise"]) - / (np.median(isovals) + results["background noise"]) - ) > ( - iqr(init_isovals) / (np.median(init_isovals) + results["background noise"]) - ): - count_initrelative += 1 + if len(init_isovals) > 0: + med_init = np.median(init_isovals) + iqr_init = iqr(init_isovals) + denom_isovals = med_isovals + results["background noise"] + denom_init = med_init + results["background noise"] + if ( + np.isfinite(denom_isovals) + and np.isfinite(denom_init) + and denom_isovals != 0 + and denom_init != 0 + ): + count_initrelative_checked += 1 + if ((iqr_isovals - results["background noise"]) / denom_isovals) > ( + iqr_init / denom_init + ): + count_initrelative += 1 + fft_denom = len(isovals) * (max(0, med_isovals) + results["background noise"]) + if not np.isfinite(fft_denom) or fft_denom <= 0: + continue f2_compare.append( np.sum(np.abs(coefs[2])) - / ( - len(isovals) - * (max(0, np.median(isovals)) + results["background noise"]) - ) + / fft_denom ) f1_compare.append( np.abs(coefs[1]) - / ( - len(isovals) - * (max(0, np.median(isovals)) + results["background noise"]) - ) + / fft_denom ) f1_compare = np.array(f1_compare) f2_compare = np.array(f2_compare) - if count_variable > (0.2 * len(checkson["R"])): + if count_skipped > 0: + logging.info( + "%s: checkfit skipped %i isophotes with too few usable samples" + % (options["ap_name"], count_skipped) + ) + if count_checked == 0: + logging.warning( + "%s: Possible failed fit! no isophotes had enough usable samples for checkfit" + % options["ap_name"] + ) + tests["isophote variability"] = False + tests["initial fit compare"] = False + tests["FFT coefficients"] = False + tests["Light symmetry"] = False + elif count_variable > (0.2 * count_checked): logging.warning( "%s: Possible failed fit! flux values highly variable along isophotes" % options["ap_name"] @@ -156,7 +195,12 @@ def Check_Fit(IMG, results, options): tests["isophote variability"] = False else: tests["isophote variability"] = True - if count_initrelative > (0.5 * len(checkson["R"])): + if count_checked == 0: + pass + elif ( + count_initrelative_checked > 0 + and count_initrelative > (0.5 * count_initrelative_checked) + ): logging.warning( "%s: Possible failed fit! flux values highly variable relative to initialization" % options["ap_name"] @@ -164,10 +208,13 @@ def Check_Fit(IMG, results, options): tests["initial fit compare"] = False else: tests["initial fit compare"] = True - if ( - np.sum(f2_compare > 0.2) > (0.1 * len(checkson["R"])) - or np.sum(f2_compare > 0.1) > (0.3 * len(checkson["R"])) - or np.sum(f2_compare > 0.05) > (0.8 * len(checkson["R"])) + if count_checked == 0: + pass + elif ( + len(f2_compare) == 0 + or np.sum(f2_compare > 0.2) > (0.1 * len(f2_compare)) + or np.sum(f2_compare > 0.1) > (0.3 * len(f2_compare)) + or np.sum(f2_compare > 0.05) > (0.8 * len(f2_compare)) ): logging.warning( "%s: Possible failed fit! poor convergence of FFT coefficients" @@ -176,10 +223,13 @@ def Check_Fit(IMG, results, options): tests["FFT coefficients"] = False else: tests["FFT coefficients"] = True - if ( - np.sum(f1_compare > 0.2) > (0.1 * len(checkson["R"])) - or np.sum(f1_compare > 0.1) > (0.3 * len(checkson["R"])) - or np.sum(f1_compare > 0.05) > (0.8 * len(checkson["R"])) + if count_checked == 0: + pass + elif ( + len(f1_compare) == 0 + or np.sum(f1_compare > 0.2) > (0.1 * len(f1_compare)) + or np.sum(f1_compare > 0.1) > (0.3 * len(f1_compare)) + or np.sum(f1_compare > 0.05) > (0.8 * len(f1_compare)) ): logging.warning( "%s: Possible failed fit! possible failed center or lopsided galaxy" @@ -190,6 +240,8 @@ def Check_Fit(IMG, results, options): tests["Light symmetry"] = True res = {"checkfit": tests} + if count_skipped > 0: + res["auxfile checkfit skipped"] = "checkfit skipped isophotes: %i" % count_skipped for t in tests: res["auxfile checkfit %s" % t] = "checkfit %s: %s" % ( t, diff --git a/src/autoprof/pipeline_steps/Isophote_Extract.py b/src/autoprof/pipeline_steps/Isophote_Extract.py index 4abe2ca9..511629e2 100644 --- a/src/autoprof/pipeline_steps/Isophote_Extract.py +++ b/src/autoprof/pipeline_steps/Isophote_Extract.py @@ -24,6 +24,7 @@ SBprof_to_COG_errorprop, _iso_extract, _iso_between, + _photutils_masked_data, LSBImage, AddLogo, _average, @@ -45,6 +46,24 @@ __all__ = ("Isophote_Extract_Forced", "Isophote_Extract", "Isophote_Extract_Photutils") +def _finite_median(values): + values = np.ma.asarray(values).compressed() + values = values[np.isfinite(values)] + if len(values) == 0: + return np.nan + return np.median(values) + + +def _finite_divide(numerator, denominator): + if ( + not np.isfinite(numerator) + or not np.isfinite(denominator) + or denominator == 0 + ): + return np.nan + return numerator / denominator + + def _Generate_Profile(IMG, results, R, parameters, options): # Create image array with background and mask applied @@ -146,6 +165,25 @@ def _Generate_Profile(IMG, results, R, parameters, options): ), sclip_nsigma=options["ap_isoclip_nsigma"] if "ap_isoclip_nsigma" in options else 5, ) + if len(isovals[0]) == 0: + pixels.append(0) + maskedpixels.append(isovals[2]) + if fluxunits == "intensity": + sb.append(np.nan) + sbE.append(np.nan) + cogdirect.append(np.nan) + else: + sb.append(99.999) + sbE.append(99.999) + cogdirect.append(99.999) + if "ap_iso_measurecoefs" in options and not options["ap_iso_measurecoefs"] is None: + measFmodes.append( + { + "a": [np.nan] * (len(options["ap_iso_measurecoefs"]) + 1), + "b": [np.nan] * (len(options["ap_iso_measurecoefs"]) + 1), + } + ) + continue isotot = np.sum(_iso_between(dat, 0, R[i], parameters[i], results["center"], mask=mask)) medflux = _average( isovals[0], @@ -996,6 +1034,7 @@ def Isophote_Extract_Photutils(IMG, results, options): SBprof_data = dict((h, []) for h in params) res = {} dat = IMG - results["background"] + photutils_dat = _photutils_masked_data(dat, results) if not "fit R" in results and not "fit photutils isolist" in results: logging.info("%s: photutils fitting and extracting image data" % options["ap_name"]) geo = EllipseGeometry( @@ -1005,13 +1044,17 @@ def Isophote_Extract_Photutils(IMG, results, options): eps=results["init ellip"], pa=results["init pa"], ) - ellipse = Photutils_Ellipse(dat, geometry=geo) + ellipse = Photutils_Ellipse(photutils_dat, geometry=geo) isolist = ellipse.fit_image(fix_center=True, linear=False) res.update( { "fit photutils isolist": isolist, - "auxfile fitlimit": "fit limit semi-major axis: %.2f pix" % isolist.sma[-1], + "auxfile fitlimit": ( + "fit limit semi-major axis: %.2f pix" % isolist.sma[-1] + if len(isolist.sma) > 0 + else "fit limit semi-major axis: no valid isophotes" + ), } ) elif not "fit photutils isolist" in results: @@ -1029,7 +1072,7 @@ def Isophote_Extract_Photutils(IMG, results, options): pa=results["fit pa"][i], ) # Extract the isophote information - ES = EllipseSample(dat, sma=results["fit R"][i], geometry=geo) + ES = EllipseSample(photutils_dat, sma=results["fit R"][i], geometry=geo) ES.update(fixed_parameters=None) list_iso.append(Isophote(ES, niter=30, valid=True, stop_code=0)) @@ -1037,7 +1080,11 @@ def Isophote_Extract_Photutils(IMG, results, options): res.update( { "fit photutils isolist": isolist, - "auxfile fitlimit": "fit limit semi-major axis: %.2f pix" % isolist.sma[-1], + "auxfile fitlimit": ( + "fit limit semi-major axis: %.2f pix" % isolist.sma[-1] + if len(isolist.sma) > 0 + else "fit limit semi-major axis: no valid isophotes" + ), } ) else: @@ -1045,27 +1092,35 @@ def Isophote_Extract_Photutils(IMG, results, options): for i in range(len(isolist.sma)): SBprof_data["R"].append(isolist.sma[i] * options["ap_pixscale"]) + medflux = _finite_median(isolist.sample[i].values[2]) if fluxunits == "intensity": SBprof_data["I"].append( - np.median(isolist.sample[i].values[2]) / options["ap_pixscale"] ** 2 + medflux / options["ap_pixscale"] ** 2 ) SBprof_data["I_e"].append(isolist.int_err[i]) SBprof_data["totflux"].append(isolist.tflux_e[i]) - SBprof_data["totflux_e"].append(isolist.rms[i] / np.sqrt(isolist.npix_e[i])) + SBprof_data["totflux_e"].append( + _finite_divide(isolist.rms[i], np.sqrt(isolist.npix_e[i])) + ) else: SBprof_data["SB"].append( - flux_to_sb( - np.median(isolist.sample[i].values[2]), - options["ap_pixscale"], - zeropoint, - ) + flux_to_sb(medflux, options["ap_pixscale"], zeropoint) + if medflux > 0 + else np.nan + ) + SBprof_data["SB_e"].append( + _finite_divide(2.5 * isolist.int_err[i], isolist.intens[i] * np.log(10)) + ) + SBprof_data["totmag"].append( + flux_to_mag(isolist.tflux_e[i], zeropoint) + if np.isfinite(isolist.tflux_e[i]) and isolist.tflux_e[i] > 0 + else np.nan ) - SBprof_data["SB_e"].append(2.5 * isolist.int_err[i] / (isolist.intens[i] * np.log(10))) - SBprof_data["totmag"].append(flux_to_mag(isolist.tflux_e[i], zeropoint)) SBprof_data["totmag_e"].append( - 2.5 - * isolist.rms[i] - / (np.sqrt(isolist.npix_e[i]) * isolist.tflux_e[i] * np.log(10)) + _finite_divide( + 2.5 * isolist.rms[i], + np.sqrt(isolist.npix_e[i]) * isolist.tflux_e[i] * np.log(10), + ) ) SBprof_data["ellip"].append(isolist.eps[i]) SBprof_data["ellip_e"].append(isolist.ellip_err[i]) @@ -1084,7 +1139,7 @@ def Isophote_Extract_Photutils(IMG, results, options): SBprof_data[k][-1] = 99.999 res.update({"prof header": params, "prof units": SBprof_units, "prof data": SBprof_data}) - if "ap_doplot" in options and options["ap_doplot"]: + if "ap_doplot" in options and options["ap_doplot"] and len(SBprof_data["R"]) > 0: if fluxunits == "intensity": Plot_I_Profile( dat, diff --git a/src/autoprof/pipeline_steps/Isophote_Fit.py b/src/autoprof/pipeline_steps/Isophote_Fit.py index 174d37b6..9ef66228 100644 --- a/src/autoprof/pipeline_steps/Isophote_Fit.py +++ b/src/autoprof/pipeline_steps/Isophote_Fit.py @@ -14,6 +14,7 @@ from ..autoprofutils.SharedFunctions import ( _iso_extract, + _photutils_masked_data, _x_to_pa, _x_to_eps, _inv_x_to_eps, @@ -37,6 +38,202 @@ ) +def _finite_isophote_samples(isovals): + isovals = np.asarray(isovals) + return isovals[np.isfinite(isovals)] + + +def _best_finite_index(losses): + losses = np.asarray(losses, dtype=float) + if not np.any(np.isfinite(losses)): + return 0 + return int(np.argmin(np.where(np.isfinite(losses), losses, np.inf))) + + +def _max_fit_mode(fit_coefs=None): + return max(fit_coefs) if not fit_coefs is None and len(fit_coefs) > 0 else 2 + + +def _has_enough_isophote_samples( + dat, + radius, + params, + center, + mask, + max_mode=2, + interp_method=None, +): + kwargs = {} + if not interp_method is None: + kwargs["interp_mask"] = False if mask is None else True + kwargs["interp_method"] = interp_method + isovals = _finite_isophote_samples( + _iso_extract( + dat, + radius, + params, + center, + mask=mask, + **kwargs, + ) + ) + return len(isovals) > 2 * max_mode + + +def _active_neighbors(i, active_mask, n): + if active_mask is None: + return [j for j in (i + 1, i - 1) if 0 <= j < n] + + active_indices = np.flatnonzero(active_mask) + neighbors = [] + inner = active_indices[active_indices < i] + outer = active_indices[active_indices > i] + if len(outer) > 0: + neighbors.append(outer[0]) + if len(inner) > 0: + neighbors.append(inner[-1]) + return neighbors + + +def _interpolate_inactive_parameters(sample_radii, parameters, active_mask): + active_indices = np.flatnonzero(active_mask) + inactive_indices = np.flatnonzero(np.logical_not(active_mask)) + if len(active_indices) < 2 or len(inactive_indices) == 0: + return + + active_radii = sample_radii[active_indices] + active_pa = np.array(list(parameters[i]["pa"] for i in active_indices)) + active_ellip = np.array(list(parameters[i]["ellip"] for i in active_indices)) + pa_s = np.sin(2 * active_pa) + pa_c = np.cos(2 * active_pa) + inv_ellip = _inv_x_to_eps(active_ellip) + + for i in inactive_indices: + radius = sample_radii[i] + parameters[i]["ellip"] = _x_to_eps(np.interp(radius, active_radii, inv_ellip)) + parameters[i]["pa"] = _x_to_pa( + ((np.arctan2(np.interp(radius, active_radii, pa_s), np.interp(radius, active_radii, pa_c))) % (2 * np.pi)) + / 2 + ) + if not parameters[i]["C"] is None: + active_C = np.array(list(parameters[j]["C"] for j in active_indices)) + parameters[i]["C"] = np.interp(radius, active_radii, active_C) + if not parameters[i]["m"] is None: + for m in range(len(parameters[i]["m"])): + active_Am = np.array(list(parameters[j]["Am"][m] for j in active_indices)) + active_Phim = np.array(list(parameters[j]["Phim"][m] for j in active_indices)) + mode = parameters[i]["m"][m] + parameters[i]["Am"][m] = np.interp(radius, active_radii, active_Am) + phase_s = np.interp(radius, active_radii, np.sin(mode * active_Phim)) + phase_c = np.interp(radius, active_radii, np.cos(mode * active_Phim)) + parameters[i]["Phim"][m] = ((np.arctan2(phase_s, phase_c)) % (2 * np.pi)) / mode + + +def _activate_inactive_radii( + dat, + sample_radii, + parameters, + active_mask, + center, + mask, + max_mode=2, + interp_method=None, +): + _interpolate_inactive_parameters(sample_radii, parameters, active_mask) + activated = [] + for i in np.flatnonzero(np.logical_not(active_mask)): + if _has_enough_isophote_samples( + dat, + sample_radii[i], + parameters[i], + center, + mask, + max_mode=max_mode, + interp_method=interp_method, + ): + active_mask[i] = True + activated.append(i) + return activated + + +def _activate_inactive_ellipse_radii(dat, sample_radii, ellip, pa, active_mask, center, mask): + active_indices = np.flatnonzero(active_mask) + inactive_indices = np.flatnonzero(np.logical_not(active_mask)) + if len(active_indices) < 2 or len(inactive_indices) == 0: + return [] + + active_radii = sample_radii[active_indices] + pa_s = np.sin(2 * pa[active_indices]) + pa_c = np.cos(2 * pa[active_indices]) + inv_ellip = _inv_x_to_eps(ellip[active_indices]) + activated = [] + for i in inactive_indices: + radius = sample_radii[i] + ellip[i] = _x_to_eps(np.interp(radius, active_radii, inv_ellip)) + pa[i] = _x_to_pa( + ((np.arctan2(np.interp(radius, active_radii, pa_s), np.interp(radius, active_radii, pa_c))) % (2 * np.pi)) + / 2 + ) + if _has_enough_isophote_samples( + dat, + radius, + {"ellip": ellip[i], "pa": pa[i]}, + center, + mask, + ): + active_mask[i] = True + activated.append(i) + return activated + + +def _determine_sample_radii( + dat, + IMG, + results, + options, + scale, + mask, + startR, + fit_limit_default, + average_func, + minR=0.0, +): + shrink = 0 + while shrink < 5: + sample_radii = [] + radius = startR + while radius < (max(IMG.shape) / 2): + sample_radii.append(radius) + isovals = _finite_isophote_samples( + _iso_extract( + dat, + radius, + {"ellip": results["init ellip"], "pa": results["init pa"]}, + results["center"], + more=False, + mask=mask, + ) + ) + if ( + len(isovals) > 0 + and radius > minR + and average_func(isovals) + < (options["ap_fit_limit"] if "ap_fit_limit" in options else fit_limit_default) + * results["background noise"] + ): + break + radius *= 1.0 + scale / (1.0 + shrink) + if len(sample_radii) < 15: + shrink += 1 + else: + break + if shrink >= 5: + raise Exception( + "Unable to initialize ellipse fit, check diagnostic plots. Possible missed center." + ) + return np.array(sample_radii) + + def Photutils_Fit(IMG, results, options): """Photutils elliptical isophote wrapper. @@ -80,7 +277,8 @@ def Photutils_Fit(IMG, results, options): eps=results["init ellip"], pa=results["init pa"], ) - ellipse = Photutils_Ellipse(dat, geometry=geo) + photutils_dat = _photutils_masked_data(dat, results) + ellipse = Photutils_Ellipse(photutils_dat, geometry=geo) isolist = ellipse.fit_image(fix_center=True, linear=False) res = { @@ -90,10 +288,14 @@ def Photutils_Fit(IMG, results, options): "fit pa": isolist.pa[1:], "fit pa_err": isolist.pa_err[1:], "fit photutils isolist": isolist, - "auxfile fitlimit": "fit limit semi-major axis: %.2f pix" % isolist.sma[-1], + "auxfile fitlimit": ( + "fit limit semi-major axis: %.2f pix" % isolist.sma[-1] + if len(isolist.sma) > 0 + else "fit limit semi-major axis: no valid isophotes" + ), } - if "ap_doplot" in options and options["ap_doplot"]: + if "ap_doplot" in options and options["ap_doplot"] and len(res["fit R"]) > 0: Plot_Isophote_Fit( dat, res["fit R"], @@ -164,33 +366,17 @@ def Isophote_Fit_FixedPhase(IMG, results, options): # Determine sampling radii ###################################################################### - shrink = 0 - while shrink < 5: - sample_radii = [max(1.0, results["psf fwhm"] / 2)] - while sample_radii[-1] < (max(IMG.shape) / 2): - isovals = _iso_extract( - dat, - sample_radii[-1], - {"ellip": results["init ellip"], "pa": results["init pa"]}, - results["center"], - more=False, - mask=mask, - ) - if ( - np.median(isovals) - < (options["ap_fit_limit"] if "ap_fit_limit" in options else 2) - * results["background noise"] - ): - break - sample_radii.append(sample_radii[-1] * (1.0 + scale / (1.0 + shrink))) - if len(sample_radii) < 15: - shrink += 1 - else: - break - if shrink >= 5: - raise Exception( - "Unable to initialize ellipse fit, check diagnostic plots. Possible missed center." - ) + sample_radii = _determine_sample_radii( + dat, + IMG, + results, + options, + scale, + mask, + max(1.0, results["psf fwhm"] / 2), + 2, + np.median, + ) ellip = np.ones(len(sample_radii)) * results["init ellip"] pa = np.ones(len(sample_radii)) * results["init pa"] logging.debug("%s: sample radii: %s" % (options["ap_name"], str(sample_radii))) @@ -230,8 +416,21 @@ def _pa_smooth(R, PA, deg): def _FFT_Robust_loss( - dat, R, PARAMS, i, C, noise, mask=None, reg_scale=1.0, robust_clip=0.15, fit_coefs=None, name="" + dat, + R, + PARAMS, + i, + C, + noise, + mask=None, + reg_scale=1.0, + robust_clip=0.15, + fit_coefs=None, + name="", + active_mask=None, ): + if fit_coefs is not None and len(fit_coefs) == 0: + fit_coefs = None isovals = _iso_extract( dat, @@ -242,62 +441,44 @@ def _FFT_Robust_loss( interp_mask=False if mask is None else True, interp_method="bicubic", ) - - try: - coefs = fft(np.clip(isovals, a_max=np.quantile(isovals, 1.0 - robust_clip), a_min=None)) - except: - coefs = np.zeros(100) - isovals = np.zeros(100) - - if fit_coefs is None: - f2_loss = np.abs(coefs[2]) / ( - len(isovals) * (max(0, np.median(isovals)) + noise / np.sqrt(len(isovals))) - ) - else: - f2_loss = np.sum(np.abs(coefs[np.array(fit_coefs)])) / ( - len(fit_coefs) - * len(isovals) - * (max(0, np.median(isovals)) + noise / np.sqrt(len(isovals))) - ) + isovals = _finite_isophote_samples(isovals) + max_mode = _max_fit_mode(fit_coefs) + f2_loss = None + if len(isovals) > 2 * max_mode: + try: + coefs = fft(np.clip(isovals, a_max=np.quantile(isovals, 1.0 - robust_clip), a_min=None)) + norm = len(isovals) * (max(0, np.median(isovals)) + noise / np.sqrt(len(isovals))) + if fit_coefs is None: + f2_loss = np.abs(coefs[2]) / norm + else: + f2_loss = np.sum(np.abs(coefs[np.array(fit_coefs)])) / (len(fit_coefs) * norm) + except Exception: + f2_loss = None reg_loss = 0 - if not PARAMS[i]["m"] is None: + has_fmodes = not PARAMS[i]["m"] is None and len(PARAMS[i]["m"]) > 0 + if has_fmodes: fmode_scale = 1.0 / len(PARAMS[i]["m"]) - if i < (len(R) - 1): - reg_loss += abs( - (PARAMS[i]["ellip"] - PARAMS[i + 1]["ellip"]) / (1 - PARAMS[i + 1]["ellip"]) - ) - reg_loss += abs(Angle_TwoAngles_sin(PARAMS[i]["pa"], PARAMS[i + 1]["pa"]) / (0.2)) - if not PARAMS[i]["m"] is None: - for m in range(len(PARAMS[i]["m"])): - reg_loss += fmode_scale * abs((PARAMS[i]["Am"][m] - PARAMS[i + 1]["Am"][m]) / 0.2) - reg_loss += fmode_scale * abs( - Angle_TwoAngles_cos( - PARAMS[i]["m"][m] * PARAMS[i]["Phim"][m], - PARAMS[i + 1]["m"][m] * PARAMS[i + 1]["Phim"][m], - ) - / (PARAMS[i]["m"][m] * 0.1) - ) - if not PARAMS[i]["C"] is None: - reg_loss += abs(np.log10(PARAMS[i]["C"] / PARAMS[i + 1]["C"])) / 0.1 - if i > 0: + for j in _active_neighbors(i, active_mask, len(R)): reg_loss += abs( - (PARAMS[i]["ellip"] - PARAMS[i - 1]["ellip"]) / (1 - PARAMS[i - 1]["ellip"]) + (PARAMS[i]["ellip"] - PARAMS[j]["ellip"]) / (1 - PARAMS[j]["ellip"]) ) - reg_loss += abs(Angle_TwoAngles_sin(PARAMS[i]["pa"], PARAMS[i - 1]["pa"]) / (0.2)) - if not PARAMS[i]["m"] is None: + reg_loss += abs(Angle_TwoAngles_sin(PARAMS[i]["pa"], PARAMS[j]["pa"]) / (0.2)) + if has_fmodes and not PARAMS[j]["m"] is None: for m in range(len(PARAMS[i]["m"])): - reg_loss += fmode_scale * abs((PARAMS[i]["Am"][m] - PARAMS[i - 1]["Am"][m]) / 0.2) + reg_loss += fmode_scale * abs((PARAMS[i]["Am"][m] - PARAMS[j]["Am"][m]) / 0.2) reg_loss += fmode_scale * abs( Angle_TwoAngles_cos( PARAMS[i]["m"][m] * PARAMS[i]["Phim"][m], - PARAMS[i - 1]["m"][m] * PARAMS[i - 1]["Phim"][m], + PARAMS[j]["m"][m] * PARAMS[j]["Phim"][m], ) / (PARAMS[i]["m"][m] * 0.1) ) - if not PARAMS[i]["C"] is None: - reg_loss += abs(np.log10(PARAMS[i]["C"] / PARAMS[i - 1]["C"])) / 0.1 + if not PARAMS[i]["C"] is None and not PARAMS[j]["C"] is None: + reg_loss += abs(np.log10(PARAMS[i]["C"] / PARAMS[j]["C"])) / 0.1 + if f2_loss is None or not np.isfinite(f2_loss): + return np.inf return f2_loss * (1 + reg_loss * reg_scale) @@ -309,44 +490,69 @@ def _FFT_Robust_Errors( E_err = np.zeros(len(R)) for ri in range(len(R)): temp_fits = [] + x0 = np.array([PARAMS[ri]["ellip"], PARAMS[ri]["pa"]], dtype=float) + if not np.all(np.isfinite(x0)): + continue for i in range(10): low_ri = max(0, ri - 1) high_ri = min(len(R) - 1, ri + 1) - temp_fits.append( - minimize( - lambda x: _FFT_Robust_loss( - dat, - [R[low_ri], R[ri] * (1 - 0.05 + i * 0.1 / 9), R[high_ri]], - [ - PARAMS[low_ri], - { - "ellip": np.clip(x[0], 0, 0.999), - "pa": x[1] % np.pi, - "m": PARAMS[ri]["m"], - "Am": PARAMS[ri]["Am"], - "Phim": PARAMS[ri]["Phim"], - }, - PARAMS[high_ri], - ], - 1, - C, - noise, - mask=mask, - reg_scale=reg_scale, - robust_clip=robust_clip, - fit_coefs=fit_coefs, - name=name, - ), - x0=[PARAMS[ri]["ellip"], PARAMS[ri]["pa"]], + trial_R = [R[low_ri], R[ri] * (1 - 0.05 + i * 0.1 / 9), R[high_ri]] + + def raw_loss(x): + if not np.all(np.isfinite(x)): + return np.inf + trial_params = deepcopy(PARAMS[ri]) + trial_params["ellip"] = float(x[0]) + trial_params["pa"] = float(x[1]) + return _FFT_Robust_loss( + dat, + trial_R, + [ + PARAMS[low_ri], + trial_params, + PARAMS[high_ri], + ], + 1, + C, + noise, + mask=mask, + reg_scale=reg_scale, + robust_clip=robust_clip, + fit_coefs=fit_coefs, + name=name, + ) + + if not np.isfinite(raw_loss(x0)): + continue + + # SLSQP expects a finite scalar objective even when a shifted isophote is unusable. + def finite_loss(x): + loss = raw_loss(x) + return float(loss) if np.isfinite(loss) else 1e30 + + try: + fit = minimize( + finite_loss, + x0=x0, method="SLSQP", + bounds=[(0, 0.999), (0, np.pi)], options={"ftol": 0.001}, - ).x + ) + except Exception: + continue + if np.all(np.isfinite(fit.x)) and np.isfinite(raw_loss(fit.x)): + temp_fits.append(fit.x) + if len(temp_fits) >= 2: + temp_fits = np.array(temp_fits) + E_err[ri] = iqr(temp_fits[:, 0], rng=[16, 84]) / 2 + PA_err[ri] = ( + Angle_Scatter(2 * (temp_fits[:, 1] % np.pi)) / 4.0 + ) # multiply by 2 to get [0, 2pi] range + elif len(temp_fits) == 1: + E_err[ri] = abs(float(temp_fits[0][0]) - PARAMS[ri]["ellip"]) + PA_err[ri] = ( + abs(Angle_TwoAngles_sin(2 * temp_fits[0][1], 2 * PARAMS[ri]["pa"])) / 2 ) - temp_fits = np.array(temp_fits) - E_err[ri] = iqr(np.clip(temp_fits[:, 0], 0, 1), rng=[16, 84]) / 2 - PA_err[ri] = ( - Angle_Scatter(2 * (temp_fits[:, 1] % np.pi)) / 4.0 - ) # multiply by 2 to get [0, 2pi] range return E_err, PA_err @@ -514,34 +720,18 @@ def Isophote_Fit_FFT_Robust(IMG, results, options): minR = 0.0 # Determine sampling radii ###################################################################### - shrink = 0 - while shrink < 5: - sample_radii = [max(1.0, results["psf fwhm"] / 2)] - while sample_radii[-1] < (max(IMG.shape) / 2): - isovals = _iso_extract( - dat, - sample_radii[-1], - {"ellip": results["init ellip"], "pa": results["init pa"]}, - results["center"], - more=False, - mask=mask, - ) - if ( - sample_radii[-1] > minR - and np.median(isovals) - < (options["ap_fit_limit"] if "ap_fit_limit" in options else 2) - * results["background noise"] - ): - break - sample_radii.append(sample_radii[-1] * (1.0 + scale / (1.0 + shrink))) - if len(sample_radii) < 15: - shrink += 1 - else: - break - if shrink >= 5: - raise Exception( - "Unable to initialize ellipse fit, check diagnostic plots. Possible missed center." - ) + sample_radii = _determine_sample_radii( + dat, + IMG, + results, + options, + scale, + mask, + max(1.0, results["psf fwhm"] / 2), + 2, + np.median, + minR=minR, + ) ellip = np.ones(len(sample_radii)) * results["init ellip"] pa = np.ones(len(sample_radii)) * results["init pa"] logging.debug("%s: sample radii: %s" % (options["ap_name"], str(sample_radii))) @@ -579,10 +769,30 @@ def Isophote_Fit_FFT_Robust(IMG, results, options): ) count_nochange = 0 use_center = copy(results["center"]) - I = np.array(range(len(sample_radii))) + active_mask = np.array( + list( + _has_enough_isophote_samples( + dat, + sample_radii[i], + parameters[i], + use_center, + mask, + max_mode=_max_fit_mode(fit_coefs), + interp_method="bicubic", + ) + for i in range(len(sample_radii)) + ) + ) + if np.sum(active_mask) < 2: + raise Exception( + "Unable to initialize ellipse fit, too few radii have enough unmasked samples." + ) param_cycle = 2 base_params = 2 + int(fit_superellipse) while count < iterlimitmax: + I = np.flatnonzero(active_mask) + if len(I) == 0: + break # Periodically include logging message if count % 10 == 0: logging.debug("%s: count: %i" % (options["ap_name"], count)) @@ -590,6 +800,7 @@ def Isophote_Fit_FFT_Robust(IMG, results, options): np.random.shuffle(I) N_perturb = int(1 + (10 / np.sqrt(count))) + activated_this_pass = False for i in I: perturbations = [] @@ -606,6 +817,7 @@ def Isophote_Fit_FFT_Robust(IMG, results, options): robust_clip=robust_clip, fit_coefs=fit_coefs, name=options["ap_name"], + active_mask=active_mask, ) for n in range(N_perturb): perturbations.append(deepcopy(parameters)) @@ -653,19 +865,39 @@ def Isophote_Fit_FFT_Robust(IMG, results, options): robust_clip=robust_clip, fit_coefs=fit_coefs, name=options["ap_name"], + active_mask=active_mask, ) - best = np.argmin(list(p[i]["loss"] for p in perturbations)) + if not np.isfinite(perturbations[0][i]["loss"]): + active_mask[i] = False + continue + losses = list(p[i]["loss"] for p in perturbations) + best = _best_finite_index(losses) if best > 0: parameters = deepcopy(perturbations[best]) del parameters[i]["loss"] count_nochange = 0 else: count_nochange += 1 + stop_threshold = iterstopnochange * max(1, len(I) - 1) if not ( - count_nochange < (iterstopnochange * (len(sample_radii) - 1)) + count_nochange < stop_threshold or count < iterlimitmin ): + activated = _activate_inactive_radii( + dat, + sample_radii, + parameters, + active_mask, + use_center, + mask, + max_mode=_max_fit_mode(fit_coefs), + interp_method="bicubic", + ) + if len(activated) > 0: + count_nochange = 0 + activated_this_pass = True + break if param_cycle > 2 or (parameters[i]["m"] is None and not fit_superellipse): break elif parameters[i]["m"] is None and fit_superellipse: @@ -678,6 +910,20 @@ def Isophote_Fit_FFT_Robust(IMG, results, options): count = 0 if fit_coefs is None: fit_coefs = (2, 4) + active_mask = np.array( + list( + _has_enough_isophote_samples( + dat, + sample_radii[ii], + parameters[ii], + use_center, + mask, + max_mode=_max_fit_mode(fit_coefs), + interp_method="bicubic", + ) + for ii in range(len(sample_radii)) + ) + ) else: logging.info( "%s: Started Fmode fitting at iteration %i" % (options["ap_name"], count) @@ -699,11 +945,27 @@ def Isophote_Fit_FFT_Robust(IMG, results, options): or options["ap_isofit_losscoefs"] is None ): fit_coefs = tuple(sorted(set([4] + list(fit_coefs)))) + active_mask = np.array( + list( + _has_enough_isophote_samples( + dat, + sample_radii[ii], + parameters[ii], + use_center, + mask, + max_mode=_max_fit_mode(fit_coefs), + interp_method="bicubic", + ) + for ii in range(len(sample_radii)) + ) + ) if ( "ap_isofit_fitcoefs_FFTinit" in options and options["ap_isofit_fitcoefs_FFTinit"] ): - for ii in I: + for ii in np.flatnonzero(active_mask): + if parameters[ii]["m"] is None or len(parameters[ii]["m"]) == 0: + continue isovals = _iso_extract( dat, sample_radii[ii], @@ -713,23 +975,20 @@ def Isophote_Fit_FFT_Robust(IMG, results, options): interp_mask=False if mask is None else True, interp_method="bicubic", ) - - if mask is None: - coefs = fft( - np.clip( - isovals, - a_max=np.quantile(isovals, 0.85), - a_min=None, - ) - ) - else: - coefs = fft( - np.clip( - isovals, - a_max=np.quantile(isovals, 0.9), - a_min=None, - ) + isovals = _finite_isophote_samples(isovals) + if len(isovals) <= 2 * max(parameters[ii]["m"]): + continue + + quantile = 0.85 if mask is None else 0.9 + coefs = fft( + np.clip( + isovals, + a_max=np.quantile(isovals, quantile), + a_min=None, ) + ) + if coefs[0] == 0 or not np.isfinite(coefs[0]): + continue for m in range(len(parameters[ii]["m"])): parameters[ii]["Am"][m] = np.abs( coefs[parameters[ii]["m"][m]] / coefs[0] @@ -738,11 +997,26 @@ def Isophote_Fit_FFT_Robust(IMG, results, options): coefs[parameters[ii]["m"][m]] ) % (2 * np.pi) - if not ( - count_nochange < (iterstopnochange * (len(sample_radii) - 1)) or count < iterlimitmin - ): + if activated_this_pass: + continue + stop_threshold = iterstopnochange * max(1, len(np.flatnonzero(active_mask)) - 1) + if not (count_nochange < stop_threshold or count < iterlimitmin): break + active_indices = np.flatnonzero(active_mask) + if len(active_indices) < 4: + raise Exception( + "Unable to fit isophotes, too few radii have enough unmasked samples." + ) + skipped_fit_radii = len(sample_radii) - len(active_indices) + if skipped_fit_radii > 0: + logging.info( + "%s: Excluded %i fit radii with too few unmasked samples" + % (options["ap_name"], skipped_fit_radii) + ) + sample_radii = np.asarray(sample_radii)[active_indices] + parameters = list(parameters[i] for i in active_indices) + logging.info("%s: Completed isohpote fit in %i itterations" % (options["ap_name"], count)) # Compute errors ###################################################################### @@ -758,7 +1032,7 @@ def Isophote_Fit_FFT_Robust(IMG, results, options): fit_coefs=fit_coefs, name=options["ap_name"], ) - for i in range(len(ellip)): + for i in range(len(parameters)): parameters[i]["ellip err"] = ellip_err[i] parameters[i]["pa err"] = pa_err[i] # Plot fitting results @@ -774,6 +1048,8 @@ def Isophote_Fit_FFT_Robust(IMG, results, options): "fit pa_err": pa_err, "auxfile fitlimit": "fit limit semi-major axis: %.2f pix" % sample_radii[-1], } + if skipped_fit_radii > 0: + res["auxfile fit masked radii"] = "fit radii excluded for insufficient samples: %i" % skipped_fit_radii if not fit_params is None: res.update({"fit Fmodes": fit_params}) for m in range(len(fit_params)): @@ -868,7 +1144,19 @@ def Isophote_Fit_Forced(IMG, results, options): ###################################################################### -def _FFT_mean_loss(dat, R, E, PA, i, C, noise, mask=None, reg_scale=1.0, name=""): +def _FFT_mean_loss( + dat, + R, + E, + PA, + i, + C, + noise, + mask=None, + reg_scale=1.0, + name="", + active_mask=None, +): isovals = _iso_extract( dat, @@ -878,25 +1166,22 @@ def _FFT_mean_loss(dat, R, E, PA, i, C, noise, mask=None, reg_scale=1.0, name="" mask=mask, interp_mask=False if mask is None else True, ) + isovals = _finite_isophote_samples(isovals) - if not np.all(np.isfinite(isovals)): - logging.warning( - "Failed to evaluate isophotal flux values, skipping this ellip/pa combination" - ) - return np.inf - - coefs = fft(isovals) - - f2_loss = np.abs(coefs[2]) / (len(isovals) * (max(0, np.mean(isovals)) + noise)) + f2_loss = None + if len(isovals) > 4: + coefs = fft(isovals) + norm = len(isovals) * (max(0, np.mean(isovals)) + noise) + if norm > 0 and np.isfinite(norm): + f2_loss = np.abs(coefs[2]) / norm reg_loss = 0 - if i < (len(R) - 1): - reg_loss += abs((E[i] - E[i + 1]) / (1 - E[i + 1])) - reg_loss += abs(Angle_TwoAngles_sin(PA[i], PA[i + 1]) / (0.3)) - if i > 0: - reg_loss += abs((E[i] - E[i - 1]) / (1 - E[i - 1])) - reg_loss += abs(Angle_TwoAngles_sin(PA[i], PA[i - 1]) / (0.3)) + for j in _active_neighbors(i, active_mask, len(R)): + reg_loss += abs((E[i] - E[j]) / (1 - E[j])) + reg_loss += abs(Angle_TwoAngles_sin(PA[i], PA[j]) / (0.3)) + if f2_loss is None or not np.isfinite(f2_loss): + return np.inf return f2_loss * (1 + reg_loss * reg_scale) @@ -965,33 +1250,17 @@ def Isophote_Fit_FFT_mean(IMG, results, options): # Determine sampling radii ###################################################################### - shrink = 0 - while shrink < 5: - sample_radii = [3 * results["psf fwhm"] / 2] - while sample_radii[-1] < (max(IMG.shape) / 2): - isovals = _iso_extract( - dat, - sample_radii[-1], - {"ellip": results["init ellip"], "pa": results["init pa"]}, - results["center"], - more=False, - mask=mask, - ) - if ( - np.mean(isovals) - < (options["ap_fit_limit"] if "ap_fit_limit" in options else 1) - * results["background noise"] - ): - break - sample_radii.append(sample_radii[-1] * (1.0 + scale / (1.0 + shrink))) - if len(sample_radii) < 15: - shrink += 1 - else: - break - if shrink >= 5: - raise Exception( - "Unable to initialize ellipse fit, check diagnostic plots. Possible missed center." - ) + sample_radii = _determine_sample_radii( + dat, + IMG, + results, + options, + scale, + mask, + 3 * results["psf fwhm"] / 2, + 1, + np.mean, + ) ellip = np.ones(len(sample_radii)) * results["init ellip"] pa = np.ones(len(sample_radii)) * results["init pa"] logging.debug("%s: sample radii: %s" % (options["ap_name"], str(sample_radii))) @@ -1006,14 +1275,34 @@ def Isophote_Fit_FFT_mean(IMG, results, options): count_nochange = 0 use_center = copy(results["center"]) - I = np.array(range(len(sample_radii))) + active_mask = np.array( + list( + _has_enough_isophote_samples( + dat, + sample_radii[i], + {"ellip": ellip[i], "pa": pa[i]}, + use_center, + mask, + ) + for i in range(len(sample_radii)) + ) + ) + if np.sum(active_mask) < 2: + raise Exception( + "Unable to initialize ellipse fit, too few radii have enough unmasked samples." + ) while count < 300 and count_nochange < (3 * len(sample_radii)): + I = np.flatnonzero(active_mask) + if len(I) == 0: + break # Periodically include logging message if count % 10 == 0: logging.debug("%s: count: %i" % (options["ap_name"], count)) count += 1 np.random.shuffle(I) + activated_this_pass = False + stop_this_pass = False for i in I: perturbations = [] perturbations.append({"ellip": copy(ellip), "pa": copy(pa)}) @@ -1028,6 +1317,7 @@ def Isophote_Fit_FFT_mean(IMG, results, options): mask=mask, reg_scale=regularize_scale if count > 4 else 0, name=options["ap_name"], + active_mask=active_mask, ) for n in range(N_perturb): perturbations.append({"ellip": copy(ellip), "pa": copy(pa)}) @@ -1051,17 +1341,57 @@ def Isophote_Fit_FFT_mean(IMG, results, options): mask=mask, reg_scale=regularize_scale if count > 4 else 0, name=options["ap_name"], + active_mask=active_mask, ) - best = np.argmin(list(p["loss"] for p in perturbations)) + if not np.isfinite(perturbations[0]["loss"]): + active_mask[i] = False + continue + losses = list(p["loss"] for p in perturbations) + best = _best_finite_index(losses) if best > 0: ellip = copy(perturbations[best]["ellip"]) pa = copy(perturbations[best]["pa"]) count_nochange = 0 else: count_nochange += 1 + stop_threshold = 3 * max(1, len(I) - 1) + if count_nochange >= stop_threshold: + activated = _activate_inactive_ellipse_radii( + dat, + sample_radii, + ellip, + pa, + active_mask, + use_center, + mask, + ) + if len(activated) > 0: + count_nochange = 0 + activated_this_pass = True + break + stop_this_pass = True + break + if activated_this_pass: + continue + if stop_this_pass: + break logging.info("%s: Completed isohpote fit in %i itterations" % (options["ap_name"], count)) + active_indices = np.flatnonzero(active_mask) + if len(active_indices) < 6: + raise Exception( + "Unable to fit isophotes, too few radii have enough unmasked samples." + ) + skipped_fit_radii = len(sample_radii) - len(active_indices) + if skipped_fit_radii > 0: + logging.info( + "%s: Excluded %i fit radii with too few unmasked samples" + % (options["ap_name"], skipped_fit_radii) + ) + sample_radii = np.asarray(sample_radii)[active_indices] + ellip = ellip[active_indices] + pa = pa[active_indices] # detect collapsed center ###################################################################### for i in range(5): @@ -1166,4 +1496,6 @@ def Isophote_Fit_FFT_mean(IMG, results, options): "fit pa_err": pa_err, "auxfile fitlimit": "fit limit semi-major axis: %.2f pix" % sample_radii[-1], } + if skipped_fit_radii > 0: + res["auxfile fit masked radii"] = "fit radii excluded for insufficient samples: %i" % skipped_fit_radii return IMG, res diff --git a/src/autoprof/pipeline_steps/Isophote_Initialize.py b/src/autoprof/pipeline_steps/Isophote_Initialize.py index 8229f855..7b1f7c0c 100644 --- a/src/autoprof/pipeline_steps/Isophote_Initialize.py +++ b/src/autoprof/pipeline_steps/Isophote_Initialize.py @@ -35,6 +35,21 @@ __all__ = ("Isophote_Init_Forced", "Isophote_Initialize", "Isophote_Initialize_mean") +def _isophote_flux(isovals): + return isovals[0] if isinstance(isovals, tuple) else isovals + + +def _has_fft_sample(isovals): + return len(_isophote_flux(isovals)) >= 3 + + +def _finite_average(values): + finite_values = [v for v in values if np.isfinite(v)] + if len(finite_values) == 0: + return np.inf + return np.mean(finite_values) + + def Isophote_Init_Forced(IMG, results, options): """Read global elliptical isophote to a galaxy from an aux file. @@ -126,10 +141,14 @@ def _fitEllip_loss(e, dat, r, p, c, n, m): mask=m, interp_mask=True, ) + if len(isovals) < 3: + return np.inf coefs = fft(np.clip(isovals, a_max=np.quantile(isovals, 0.85), a_min=None)) - return (iqr(isovals, rng=[16, 84]) / 2 + np.abs(coefs[2]) / len(isovals)) / ( - max(0, np.median(isovals)) + n - ) + denominator = max(0, np.median(isovals)) + n + if (not np.isfinite(denominator)) or denominator <= 0: + return np.inf + loss = (iqr(isovals, rng=[16, 84]) / 2 + np.abs(coefs[2]) / len(isovals)) / denominator + return loss if np.isfinite(loss) else np.inf def Isophote_Initialize(IMG, results, options): @@ -201,12 +220,12 @@ def Isophote_Initialize(IMG, results, options): mask = None if "ap_isoinit_R_set" in options: - circ_ellipse_radii = np.logspace( + sample_radii = np.logspace( np.log10(circ_ellipse_radii[0]), np.log10(options["ap_isoinit_R_set"]), 10, - ).tolist() - for r in circ_ellipse_radii[1:]: + ) + for r in sample_radii[1:]: isovals = _iso_extract( dat, r, @@ -218,14 +237,18 @@ def Isophote_Initialize(IMG, results, options): sclip_nsigma=3, interp_mask=True, ) + if not _has_fft_sample(isovals): + continue + circ_ellipse_radii.append(r) coefs = fft(isovals[0]) allphase.append(coefs[2]) else: - while circ_ellipse_radii[-1] < (len(IMG) / 2): - circ_ellipse_radii.append(circ_ellipse_radii[-1] * (1 + 0.2)) + r = circ_ellipse_radii[-1] + while r < (len(IMG) / 2): + r *= 1 + 0.2 isovals = _iso_extract( dat, - circ_ellipse_radii[-1], + r, {"ellip": 0.0, "pa": 0.0}, results["center"], more=True, @@ -234,6 +257,9 @@ def Isophote_Initialize(IMG, results, options): sclip_nsigma=3, interp_mask=True, ) + if not _has_fft_sample(isovals): + continue + circ_ellipse_radii.append(r) coefs = fft(isovals[0]) allphase.append(coefs[2]) # Stop when at 3 time background noise @@ -247,6 +273,9 @@ def Isophote_Initialize(IMG, results, options): ): break + if len(allphase) == 0 or len(circ_ellipse_radii) < 2: + raise ValueError("Could not initialize isophotes: no finite samples found.") + logging.info("%s: init scale: %f pix" % (options["ap_name"], circ_ellipse_radii[-1])) # Find global position angle. phase = (-Angle_Median(np.angle(allphase[-5:])) / 2) % np.pi @@ -258,49 +287,47 @@ def Isophote_Initialize(IMG, results, options): test_f2 = [] for e in test_ellip: test_f2.append( - sum( - list( - _fitEllip_loss( - e, - dat, - circ_ellipse_radii[-2] * m, - phase, - results["center"], - results["background noise"], - mask, - ) - for m in np.linspace(0.8, 1.2, 5) + _finite_average( + _fitEllip_loss( + e, + dat, + circ_ellipse_radii[-2] * m, + phase, + results["center"], + results["background noise"], + mask, ) + for m in np.linspace(0.8, 1.2, 5) ) ) + if not np.any(np.isfinite(test_f2)): + raise ValueError("Could not initialize isophotes: no finite ellipticity samples found.") # Find global ellipticity: second pass ellip = test_ellip[np.argmin(test_f2)] test_ellip = np.linspace(ellip - 0.05, ellip + 0.05, 15) test_f2 = [] for e in test_ellip: test_f2.append( - sum( - list( - _fitEllip_loss( - e, - dat, - circ_ellipse_radii[-2] * m, - phase, - results["center"], - results["background noise"], - mask, - ) - for m in np.linspace(0.8, 1.2, 5) + _finite_average( + _fitEllip_loss( + e, + dat, + circ_ellipse_radii[-2] * m, + phase, + results["center"], + results["background noise"], + mask, ) + for m in np.linspace(0.8, 1.2, 5) ) ) + if not np.any(np.isfinite(test_f2)): + raise ValueError("Could not initialize isophotes: no finite ellipticity samples found.") ellip = test_ellip[np.argmin(test_f2)] res = minimize( - lambda e, d, r, p, c, n, msk: sum( - list( + lambda e, d, r, p, c, n, msk: _finite_average( _fitEllip_loss(_x_to_eps(e[0]), d, r * m, p, c, n, msk) for m in np.linspace(0.8, 1.2, 5) - ) ), x0=_inv_x_to_eps(ellip), args=( @@ -319,7 +346,7 @@ def Isophote_Initialize(IMG, results, options): ] }, ) - if res.success: + if res.success and np.isfinite(res.fun): logging.debug( "%s: using optimal ellipticity %.3f over grid ellipticity %.3f" % (options["ap_name"], _x_to_eps(res.x[0]), ellip) @@ -336,6 +363,7 @@ def Isophote_Initialize(IMG, results, options): 10, ) errallphase = [] + err_radii = [] for rr in RR: isovals = _iso_extract( dat, @@ -343,16 +371,24 @@ def Isophote_Initialize(IMG, results, options): {"ellip": 0.0, "pa": 0.0}, results["center"], more=True, + mask=mask, sigmaclip=True, sclip_nsigma=3, interp_mask=True, ) + if not _has_fft_sample(isovals): + continue coefs = fft(isovals[0]) errallphase.append(coefs[2]) - sample_pas = (-np.angle(1j * np.array(errallphase) / np.mean(errallphase)) / 2) % np.pi - pa_err = iqr(sample_pas, rng=[16, 84]) / 2 - res_multi = map( - lambda rrp: minimize( + err_radii.append(rr) + if len(errallphase) > 0 and np.isfinite(np.mean(errallphase)) and np.mean(errallphase) != 0: + sample_pas = (-np.angle(1j * np.array(errallphase) / np.mean(errallphase)) / 2) % np.pi + pa_err = iqr(sample_pas, rng=[16, 84]) / 2 + else: + sample_pas = np.array([]) + pa_err = np.nan + res_multi = [ + minimize( lambda e, d, r, p, c, n, m: _fitEllip_loss(_x_to_eps(e[0]), d, r, p, c, n, m), x0=_inv_x_to_eps(ellip), args=( @@ -370,10 +406,15 @@ def Isophote_Initialize(IMG, results, options): [_inv_x_to_eps(ellip) + 1 / 15], ] }, - ), - zip(RR, sample_pas), - ) - ellip_err = iqr(list(_x_to_eps(rm.x[0]) for rm in res_multi), rng=[16, 84]) / 2 + ) + for rrp in zip(err_radii, sample_pas) + ] + ellip_samples = [ + _x_to_eps(rm.x[0]) + for rm in res_multi + if rm.success and np.isfinite(rm.fun) and np.all(np.isfinite(rm.x)) + ] + ellip_err = iqr(ellip_samples, rng=[16, 84]) / 2 if len(ellip_samples) > 0 else np.nan circ_ellipse_radii = np.array(circ_ellipse_radii) @@ -409,10 +450,16 @@ def Isophote_Initialize(IMG, results, options): } -def _fitEllip_mean_loss(e, dat, r, p, c, n): - isovals = _iso_extract(dat, r, {"ellip": e, "pa": p}, c) +def _fitEllip_mean_loss(e, dat, r, p, c, n, m): + isovals = _iso_extract(dat, r, {"ellip": e, "pa": p}, c, mask=m, interp_mask=True) + if len(isovals) < 3: + return np.inf coefs = fft(isovals) - return np.abs(coefs[2]) / (len(isovals) * (max(0, np.mean(isovals)) + n)) + denominator = len(isovals) * (max(0, np.mean(isovals)) + n) + if (not np.isfinite(denominator)) or denominator <= 0: + return np.inf + loss = np.abs(coefs[2]) / denominator + return loss if np.isfinite(loss) else np.inf def Isophote_Initialize_mean(IMG, results, options): @@ -460,21 +507,33 @@ def Isophote_Initialize_mean(IMG, results, options): circ_ellipse_radii = [results["psf fwhm"]] allphase = [] dat = IMG - results["background"] + mask = results["mask"] if "mask" in results else None + if not np.any(mask): + mask = None - while circ_ellipse_radii[-1] < (len(IMG) / 2): - circ_ellipse_radii.append(circ_ellipse_radii[-1] * (1 + 0.2)) + r = circ_ellipse_radii[-1] + while r < (len(IMG) / 2): + r *= 1 + 0.2 isovals = _iso_extract( dat, - circ_ellipse_radii[-1], + r, {"ellip": 0.0, "pa": 0.0}, results["center"], more=True, + mask=mask, + interp_mask=True, ) + if not _has_fft_sample(isovals): + continue + circ_ellipse_radii.append(r) coefs = fft(isovals[0]) allphase.append(coefs[2]) # Stop when at 3 times background noise if np.mean(isovals[0]) < (3 * results["background noise"]) and len(circ_ellipse_radii) > 4: break + if len(allphase) == 0 or len(circ_ellipse_radii) < 2: + raise ValueError("Could not initialize isophotes: no finite samples found.") + logging.info("%s: init scale: %f pix" % (options["ap_name"], circ_ellipse_radii[-1])) # Find global position angle. phase = ( @@ -486,47 +545,47 @@ def Isophote_Initialize_mean(IMG, results, options): test_f2 = [] for e in test_ellip: test_f2.append( - sum( - list( - _fitEllip_mean_loss( - e, - dat, - circ_ellipse_radii[-2] * m, - phase, - results["center"], - results["background noise"], - ) - for m in np.linspace(0.8, 1.2, 5) + _finite_average( + _fitEllip_mean_loss( + e, + dat, + circ_ellipse_radii[-2] * m, + phase, + results["center"], + results["background noise"], + mask, ) + for m in np.linspace(0.8, 1.2, 5) ) ) + if not np.any(np.isfinite(test_f2)): + raise ValueError("Could not initialize isophotes: no finite ellipticity samples found.") # Find global ellipticity: second pass ellip = test_ellip[np.argmin(test_f2)] test_ellip = np.linspace(ellip - 0.05, ellip + 0.05, 15) test_f2 = [] for e in test_ellip: test_f2.append( - sum( - list( - _fitEllip_mean_loss( - e, - dat, - circ_ellipse_radii[-2] * m, - phase, - results["center"], - results["background noise"], - ) - for m in np.linspace(0.8, 1.2, 5) + _finite_average( + _fitEllip_mean_loss( + e, + dat, + circ_ellipse_radii[-2] * m, + phase, + results["center"], + results["background noise"], + mask, ) + for m in np.linspace(0.8, 1.2, 5) ) ) + if not np.any(np.isfinite(test_f2)): + raise ValueError("Could not initialize isophotes: no finite ellipticity samples found.") ellip = test_ellip[np.argmin(test_f2)] res = minimize( - lambda e, d, r, p, c, n: sum( - list( - _fitEllip_mean_loss(_x_to_eps(e[0]), d, r * m, p, c, n) + lambda e, d, r, p, c, n, msk: _finite_average( + _fitEllip_mean_loss(_x_to_eps(e[0]), d, r * m, p, c, n, msk) for m in np.linspace(0.8, 1.2, 5) - ) ), x0=_inv_x_to_eps(ellip), args=( @@ -535,6 +594,7 @@ def Isophote_Initialize_mean(IMG, results, options): phase, results["center"], results["background noise"], + mask, ), method="Nelder-Mead", options={ @@ -544,7 +604,7 @@ def Isophote_Initialize_mean(IMG, results, options): ] }, ) - if res.success: + if res.success and np.isfinite(res.fun): logging.debug( "%s: using optimal ellipticity %.3f over grid ellipticity %.3f" % (options["ap_name"], _x_to_eps(res.x[0]), ellip) @@ -559,17 +619,33 @@ def Isophote_Initialize_mean(IMG, results, options): 10, ) errallphase = [] + err_radii = [] for rr in RR: - isovals = _iso_extract(dat, rr, {"ellip": 0.0, "pa": 0.0}, results["center"], more=True) + isovals = _iso_extract( + dat, + rr, + {"ellip": 0.0, "pa": 0.0}, + results["center"], + more=True, + mask=mask, + interp_mask=True, + ) + if not _has_fft_sample(isovals): + continue coefs = fft(isovals[0]) errallphase.append(coefs[2]) - sample_pas = (-np.angle(1j * np.array(errallphase) / np.mean(errallphase)) / 2) % np.pi - pa_err = np.std(sample_pas) - res_multi = map( - lambda rrp: minimize( - lambda e, d, r, p, c, n: _fitEllip_mean_loss(_x_to_eps(e[0]), d, r, p, c, n), + err_radii.append(rr) + if len(errallphase) > 0 and np.isfinite(np.mean(errallphase)) and np.mean(errallphase) != 0: + sample_pas = (-np.angle(1j * np.array(errallphase) / np.mean(errallphase)) / 2) % np.pi + pa_err = np.std(sample_pas) + else: + sample_pas = np.array([]) + pa_err = np.nan + res_multi = [ + minimize( + lambda e, d, r, p, c, n, m: _fitEllip_mean_loss(_x_to_eps(e[0]), d, r, p, c, n, m), x0=_inv_x_to_eps(ellip), - args=(dat, rrp[0], rrp[1], results["center"], results["background noise"]), + args=(dat, rrp[0], rrp[1], results["center"], results["background noise"], mask), method="Nelder-Mead", options={ "initial_simplex": [ @@ -577,10 +653,15 @@ def Isophote_Initialize_mean(IMG, results, options): [_inv_x_to_eps(ellip) + 1 / 15], ] }, - ), - zip(RR, sample_pas), - ) - ellip_err = np.std(list(_x_to_eps(rm.x[0]) for rm in res_multi)) + ) + for rrp in zip(err_radii, sample_pas) + ] + ellip_samples = [ + _x_to_eps(rm.x[0]) + for rm in res_multi + if rm.success and np.isfinite(rm.fun) and np.all(np.isfinite(rm.x)) + ] + ellip_err = np.std(ellip_samples) if len(ellip_samples) > 0 else np.nan circ_ellipse_radii = np.array(circ_ellipse_radii) diff --git a/src/autoprof/pipeline_steps/Radial_Profiles.py b/src/autoprof/pipeline_steps/Radial_Profiles.py index a82d5e22..68ffe5da 100644 --- a/src/autoprof/pipeline_steps/Radial_Profiles.py +++ b/src/autoprof/pipeline_steps/Radial_Profiles.py @@ -198,7 +198,12 @@ def Radial_Profiles(IMG, results, options): ) ) isovals[1] -= pa[i] - avgmedflux = [] + radius_medflux = [] + if len(isovals[0]) == 0: + for sa_i in range(len(wedgeangles)): + sb[sa_i].append(99.999) + sbE[sa_i].append(99.999) + continue for sa_i in range(len(wedgeangles)): aselect = np.abs(Angle_TwoAngles_cos(wedgeangles[sa_i], isovals[1])) < ( @@ -214,7 +219,7 @@ def Radial_Profiles(IMG, results, options): if "ap_isoaverage_method" in options else "median", ) - avgmedflux.append(medflux) + radius_medflux.append(medflux) scatflux = _scatter( isovals[0][aselect], options["ap_isoaverage_method"] @@ -231,7 +236,9 @@ def Radial_Profiles(IMG, results, options): if medflux > 0 else 99.999 ) - avgmedflux = np.mean(avgmedflux) + # Only update the rolling flux state from radii with valid wedge samples. + if len(radius_medflux) > 0: + avgmedflux = np.mean(radius_medflux) if "prof header" in results: newprofheader = results["prof header"] diff --git a/src/autoprof/pipeline_steps/Slice_Profiles.py b/src/autoprof/pipeline_steps/Slice_Profiles.py index 78c4c455..55995e4c 100644 --- a/src/autoprof/pipeline_steps/Slice_Profiles.py +++ b/src/autoprof/pipeline_steps/Slice_Profiles.py @@ -93,7 +93,7 @@ def Slice_Profile(IMG, results, options): """ - dat = IMG - (results["background"] if "background" in results else np.median(IMG)) + dat = IMG - (results["background"] if "background" in results else np.nanmedian(IMG)) zeropoint = options["ap_zeropoint"] if "ap_zeropoint" in options else 22.5 use_anchor = ( @@ -147,7 +147,8 @@ def Slice_Profile(IMG, results, options): % (options["ap_name"], use_step) ) - F, X = _iso_line(dat, use_length, use_width, use_pa, use_anchor, more=False) + mask = results["mask"] if "mask" in results else None + F, X = _iso_line(dat, use_length, use_width, use_pa, use_anchor, more=False, mask=mask) windows = np.arange(0, use_length, use_step) @@ -158,7 +159,16 @@ def Slice_Profile(IMG, results, options): sb_sclip_e = [] for i in range(len(windows) - 1): isovals = F[np.logical_and(X >= windows[i], X < windows[i + 1])] - isovals_sclip = Sigma_Clip_Upper(isovals, iterations=10, nsigma=5) + if len(isovals) == 0: + sb.append(99.999) + sb_e.append(99.999) + sb_sclip.append(99.999) + sb_sclip_e.append(99.999) + continue + sclim = Sigma_Clip_Upper(isovals, iterations=10, nsigma=5) + isovals_sclip = isovals[isovals < sclim] + if len(isovals_sclip) == 0: + isovals_sclip = isovals medflux = _average( isovals, @@ -204,7 +214,7 @@ def Slice_Profile(IMG, results, options): ( 2.5 * scatflux_sclip - / (np.sqrt(len(isovals)) * medflux_sclip * np.log(10)) + / (np.sqrt(len(isovals_sclip)) * medflux_sclip * np.log(10)) ) if medflux_sclip > 0 else 99.999 @@ -221,24 +231,25 @@ def Slice_Profile(IMG, results, options): f.write( "# flux sum: %f\n" % (np.sum(F[np.logical_and(X >= 0, X <= use_length)])) ) + profile_values = F[np.logical_and(X >= 0, X <= use_length)] f.write( "# flux mean: %f\n" - % (_average(F[np.logical_and(X >= 0, X <= use_length)], "mean")) + % (_average(profile_values, "mean") if len(profile_values) > 0 else np.nan) ) f.write( "# flux median: %f\n" - % (_average(F[np.logical_and(X >= 0, X <= use_length)], "median")) + % (_average(profile_values, "median") if len(profile_values) > 0 else np.nan) ) f.write( "# flux mode: %f\n" - % (_average(F[np.logical_and(X >= 0, X <= use_length)], "mode")) + % (_average(profile_values, "mode") if len(profile_values) > 0 else np.nan) ) f.write( - "# flux std: %f\n" % (np.std(F[np.logical_and(X >= 0, X <= use_length)])) + "# flux std: %f\n" % (np.std(profile_values) if len(profile_values) > 0 else np.nan) ) f.write( "# flux 16-84%% range: %f\n" - % (iqr(F[np.logical_and(X >= 0, X <= use_length)], rng=[16, 84])) + % (iqr(profile_values, rng=[16, 84]) if len(profile_values) > 0 else np.nan) ) f.write("R,sb,sb_e,sb_sclip,sb_sclip_e\n") f.write("arcsec,mag*arcsec^-2,mag*arcsec^-2,mag*arcsec^-2,mag*arcsec^-2\n")