From 731527c8775f6685ea2cf1773970992e4eb6243e Mon Sep 17 00:00:00 2001 From: karpov-sv Date: Mon, 2 Feb 2026 15:35:55 +0100 Subject: [PATCH 1/3] Add annulus stats with robust sigma clipping --- sep.pyx | 167 ++++++++++++++++++++++++- src/aperture.c | 330 +++++++++++++++++++++++++++++++++++++++++++++++++ src/sep.h | 25 ++++ test.py | 65 ++++++++++ 4 files changed, 585 insertions(+), 2 deletions(-) diff --git a/sep.pyx b/sep.pyx index d76505c..be241ec 100644 --- a/sep.pyx +++ b/sep.pyx @@ -165,6 +165,13 @@ cdef extern from "sep.h": int id, int subpix, short inflags, double *sum, double *sumerr, double *area, short *flag) + int sep_stats_circann(const sep_image *image, + double x, double y, double rin, double rout, + int id, int subpix, short inflags, + double clip_sigma, int clip_iters, + double *mean, double *std, double *median, + double *mad_std, double *mean_clip, short *flag) + int sep_sum_ellipse(const sep_image *image, double x, double y, double a, double b, double theta, double r, int id, int subpix, short inflags, @@ -952,12 +959,12 @@ def sum_circle(np.ndarray data not None, x, y, r, flags : `~numpy.ndarray` Integer giving flags. (0 if no flags set.) + """ cdef double flux1, fluxerr1, area1, cdef double bkgflux, bkgfluxerr, bkgarea cdef short flag1, bkgflag - cdef size_t i cdef int status cdef np.broadcast it cdef sep_image im @@ -1151,7 +1158,6 @@ def sum_circann(np.ndarray data not None, x, y, rin, rout, """ cdef double area1 - cdef size_t i cdef int status cdef np.broadcast it cdef sep_image im @@ -1212,6 +1218,163 @@ def sum_circann(np.ndarray data not None, x, y, rin, rout, return sum, sumerr, flag +@cython.boundscheck(False) +@cython.wraparound(False) +def stats_circann(np.ndarray data not None, x, y, rin, rout, + var=None, err=None, gain=None, np.ndarray mask=None, + double maskthresh=0.0, seg_id=None, np.ndarray segmap=None, + int subpix=5, double clip_sigma=3.0, int clip_iters=5): + """stats_circann(data, x, y, rin, rout, err=None, var=None, mask=None, + maskthresh=0.0, seg_id=None, segmap=None, gain=None, + subpix=5) + + Compute statistics in circular annular aperture(s). + + Parameters + ---------- + data : `~numpy.ndarray` + 2-d array to be sampled. + + x, y, rin, rout : array_like + Center coordinates and inner and outer radii of aperture(s). + ``x`` corresponds to the second ("fast") axis of the input array + and ``y`` corresponds to the first ("slow") axis. + ``x, y = (0.0, 0.0)`` corresponds to the center of the first + element of the array. These inputs obey numpy broadcasting rules. + It is required that ``rout >= rin >= 0.0``. + + err, var : float or ndarray + Error *or* variance (specify at most one). Accepted for API + compatibility; not used in the statistics. + + mask : `~numpy.ndarray`, optional + Mask array. If supplied, a given pixel is masked if its value + is greater than ``maskthresh``. Masked pixels are ignored (no + correction applied). + + maskthresh : float, optional + Threshold for a pixel to be masked. Default is ``0.0``. + + segmap : `~numpy.ndarray`, optional + Segmentation image with dimensions of ``data`` and dtype ``np.int32``. + This is an optional input and corresponds to the segmentation map + output by `~sep.extract`. + + seg_id : array_like, optional + Array of segmentation ids used to mask additional pixels in the image. + Dimensions correspond to the dimensions of ``x`` and ``y``. The + behavior differs depending on whether ``seg_id`` is negative or + positive. If ``seg_id`` is positive, all pixels belonging to other + objects are masked. (Pixel ``j, i`` is masked if ``seg[j, i] != seg_id + and seg[j, i] != 0``). If ``seg_id`` is negative, all pixels other + than those belonging to the object of interest are masked. (Pixel ``j, + i`` is masked if ``seg[j, i] != -seg_id``). NB: must be included if + ``segmap`` is provided. + + gain : float, optional + Accepted for API compatibility; not used in the statistics. + + subpix : int, optional + Subpixel sampling factor. Default is 5. If 0, exact overlap is used. + + clip_sigma : float, optional + Sigma value for clipping. Default is 3.0. + + clip_iters : int, optional + Maximum number of clipping iterations. Default is 5. + + Returns + ------- + mean : `~numpy.ndarray` + Weighted mean of the data within the aperture. + + std : `~numpy.ndarray` + Weighted standard deviation of the data within the aperture. + + median : `~numpy.ndarray` + Weighted median of the data within the aperture. + + mad_std : `~numpy.ndarray` + Weighted median absolute deviation scaled by 1.4826, providing a + robust estimate of the standard deviation for Gaussian data. + + mean_clip : `~numpy.ndarray` + Sigma-clipped mean using median/MAD. + + flags : `~numpy.ndarray` + Integer giving flags. (0 if no flags set.) + + Notes + ----- + If no unmasked pixels fall within the annulus, outputs are ``NaN`` and + the ``APER_ALLMASKED`` flag is set. + """ + + cdef int status + cdef np.broadcast it + cdef sep_image im + + # Test for segmap without seg_id. Nothing happens if seg_id supplied but + # without segmap. + if (segmap is not None) and (seg_id is None): + raise ValueError('`segmap` supplied but not `seg_id`.') + + _parse_arrays(data, err, var, mask, segmap, &im) + im.maskthresh = maskthresh + if gain is not None: + im.gain = gain + + # convert inputs to double arrays + dt = np.dtype(np.double) + dint = np.dtype(np.int32) + x = np.require(x, dtype=dt) + y = np.require(y, dtype=dt) + rin = np.require(rin, dtype=dt) + rout = np.require(rout, dtype=dt) + + # Segmentation image and ids with same dimensions as x, y, etc. + if seg_id is not None: + seg_id = np.require(seg_id, dtype=dint) + if seg_id.shape != x.shape: + raise ValueError('Shapes of `x` and `seg_id` do not match') + else: + seg_id = np.zeros(len(x), dtype=dint) + + # allocate output arrays + shape = np.broadcast(x, y, rin, rout).shape + mean = np.empty(shape, dt) + std = np.empty(shape, dt) + median = np.empty(shape, dt) + mad_std = np.empty(shape, dt) + mean_clip = np.empty(shape, dt) + flag = np.empty(shape, np.short) + + it = np.broadcast(x, y, rin, rout, seg_id, mean, std, median, mad_std, mean_clip, flag) + + while np.PyArray_MultiIter_NOTDONE(it): + status = sep_stats_circann( + &im, + (np.PyArray_MultiIter_DATA(it, 0))[0], + (np.PyArray_MultiIter_DATA(it, 1))[0], + (np.PyArray_MultiIter_DATA(it, 2))[0], + (np.PyArray_MultiIter_DATA(it, 3))[0], + (np.PyArray_MultiIter_DATA(it, 4))[0], + subpix, SEP_MASK_IGNORE, + clip_sigma, clip_iters, + np.PyArray_MultiIter_DATA(it, 5), + np.PyArray_MultiIter_DATA(it, 6), + np.PyArray_MultiIter_DATA(it, 7), + np.PyArray_MultiIter_DATA(it, 8), + np.PyArray_MultiIter_DATA(it, 9), + np.PyArray_MultiIter_DATA(it, 10)) + + _assert_ok(status) + + np.PyArray_MultiIter_NEXT(it) + + return mean, std, median, mad_std, mean_clip, flag + + def sum_ellipse(np.ndarray data not None, x, y, a, b, theta, r=1.0, var=None, err=None, gain=None, np.ndarray mask=None, double maskthresh=0.0, diff --git a/src/aperture.c b/src/aperture.c index ef8b8f4..b9883c6 100644 --- a/src/aperture.c +++ b/src/aperture.c @@ -184,6 +184,41 @@ static void oversamp_ann_ellipse(double r, double b, double * r_in2, double * r_ *r_out2 = (*r_out2) * (*r_out2); } +typedef struct { + double v; + double w; +} valweight; + +static int cmp_valweight(const void *a, const void *b) { + double va = ((const valweight *)a)->v; + double vb = ((const valweight *)b)->v; + if (va < vb) { + return -1; + } + if (va > vb) { + return 1; + } + return 0; +} + +static double weighted_median_sorted(const valweight *arr, int64_t n, double totw) { + double acc = 0.0; + int64_t i; + + if (n == 0 || totw <= 0.0) { + return NAN; + } + + for (i = 0; i < n; i++) { + acc += arr[i].w; + if (acc >= 0.5 * totw) { + return arr[i].v; + } + } + + return arr[n - 1].v; +} + /*****************************************************************************/ /* circular aperture */ @@ -295,6 +330,301 @@ static void oversamp_ann_ellipse(double r, double b, double * r_in2, double * r_ #undef APER_COMPARE2 #undef APER_COMPARE3 +int sep_stats_circann( + const sep_image * im, + double x, + double y, + double rin, + double rout, + int id, + int subpix, + short inflag, + double clip_sigma, + int clip_iters, + double * mean, + double * std, + double * median, + double * mad_std, + double * mean_clip, + short * flag +) { + PIXTYPE pix; + double dx, dy, dx1, dy2, offset, scale, scale2, rpix2, overlap; + double rin2, rin_in2, rin_out2, rout2, rout_in2, rout_out2; + double totw, sumw, varw, diff, wsum, vsum; + double med, madv, sig, totw_keep; + double lo, hi; + int64_t ix, iy, i, xmin, xmax, ymin, ymax, sx, sy, pos, size, msize, ssize, nkeep; + int64_t nvals, cap; + int status, ismasked, changed, iter; + const BYTE *datat, *maskt, *segt; + converter convert, mconvert, sconvert; + valweight *vw = NULL; + valweight *tmp = NULL; + char *keep = NULL; + const double mad_scale = 1.4826; + + /* input checks */ + if (!(rin >= 0.0 && rout >= rin)) { + return ILLEGAL_APER_PARAMS; + } + if (subpix < 0) { + return ILLEGAL_SUBPIX; + } + if (clip_sigma <= 0.0 || clip_iters < 0) { + return ILLEGAL_APER_PARAMS; + } + + (void)inflag; + + *flag = 0; + *mean = NAN; + *std = NAN; + *median = NAN; + *mad_std = NAN; + *mean_clip = NAN; + + rin2 = rin * rin; + rout2 = rout * rout; + oversamp_ann_circle(rin, &rin_in2, &rin_out2); + oversamp_ann_circle(rout, &rout_in2, &rout_out2); + + if (subpix > 0) { + scale = 1.0 / subpix; + scale2 = scale * scale; + offset = 0.5 * (scale - 1.0); + } else { + scale = 0.0; + scale2 = 0.0; + offset = 0.0; + } + + maskt = NULL; + segt = NULL; + msize = 0; + ssize = 0; + + /* get data converter(s) for input array(s) */ + if ((status = get_converter(im->dtype, &convert, &size))) { + return status; + } + if (im->mask && (status = get_converter(im->mdtype, &mconvert, &msize))) { + return status; + } + if (im->segmap && (status = get_converter(im->sdtype, &sconvert, &ssize))) { + return status; + } + + /* get extent of box */ + boxextent(x, y, rout, rout, im->w, im->h, &xmin, &xmax, &ymin, &ymax, flag); + + cap = (xmax - xmin) * (ymax - ymin); + if (cap <= 0) { + *flag |= SEP_APER_ALLMASKED; + return RETURN_OK; + } + + status = RETURN_OK; + QMALLOC(vw, valweight, cap, status); + + nvals = 0; + totw = 0.0; + + /* loop over rows in the box */ + for (iy = ymin; iy < ymax; iy++) { + /* set pointers to the start of this row */ + pos = (iy % im->h) * im->w + xmin; + datat = MSVC_VOID_CAST im->data + pos * size; + if (im->mask) { + maskt = MSVC_VOID_CAST im->mask + pos * msize; + } + if (im->segmap) { + segt = MSVC_VOID_CAST im->segmap + pos * ssize; + } + + /* loop over pixels in this row */ + for (ix = xmin; ix < xmax; ix++) { + dx = ix - x; + dy = iy - y; + rpix2 = dx * dx + dy * dy; + if ((rpix2 < rout_out2) && (rpix2 > rin_in2)) { + if ((rpix2 > rout_in2) || (rpix2 < rin_out2)) { + if (subpix == 0) { + overlap = circoverlap(dx - 0.5, dy - 0.5, dx + 0.5, dy + 0.5, rout) + - circoverlap(dx - 0.5, dy - 0.5, dx + 0.5, dy + 0.5, rin); + } else { + dx += offset; + dy += offset; + overlap = 0.0; + for (sy = subpix; sy--; dy += scale) { + dx1 = dx; + dy2 = dy * dy; + for (sx = subpix; sx--; dx1 += scale) { + rpix2 = dx1 * dx1 + dy2; + if ((rpix2 < rout2) && (rpix2 > rin2)) { + overlap += scale2; + } + } + } + } + } else { + overlap = 1.0; + } + + if (overlap > 0.0) { + pix = convert(datat); + + ismasked = 0; + if (im->mask && (mconvert(maskt) > im->maskthresh)) { + ismasked = 1; + } + + /* Segmentation image: + + If `id` is negative, require segmented pixels within the + aperture. + + If `id` is positive, mask pixels with nonzero segment ids + not equal to `id`. + + */ + if (im->segmap) { + if (id > 0) { + if ((sconvert(segt) > 0.) && (sconvert(segt) != id)) { + ismasked = 1; + } + } else { + if (sconvert(segt) != -1 * id) { + ismasked = 1; + } + } + } + + if (ismasked) { + *flag |= SEP_APER_HASMASKED; + } else { + vw[nvals].v = pix; + vw[nvals].w = overlap; + totw += overlap; + nvals++; + } + } + } + + /* increment pointers by one element */ + datat += size; + maskt += msize; + segt += ssize; + } + } + + if (nvals == 0 || totw <= 0.0) { + *flag |= SEP_APER_ALLMASKED; + goto exit; + } + + /* mean and standard deviation */ + sumw = 0.0; + for (i = 0; i < nvals; i++) { + if (vw[i].w > 0.0) { + sumw += vw[i].w * vw[i].v; + } + } + *mean = sumw / totw; + + varw = 0.0; + for (i = 0; i < nvals; i++) { + if (vw[i].w > 0.0) { + diff = vw[i].v - *mean; + varw += vw[i].w * diff * diff; + } + } + *std = sqrt(varw / totw); + + /* sigma-clipped mean using median/MAD (3-sigma, max 5 iterations) */ + status = RETURN_OK; + QMALLOC(keep, char, nvals, status); + QMALLOC(tmp, valweight, nvals, status); + for (i = 0; i < nvals; i++) { + keep[i] = 1; + } + + changed = 1; + for (iter = 0; iter < clip_iters; iter++) { + nkeep = 0; + totw_keep = 0.0; + for (i = 0; i < nvals; i++) { + if (keep[i] && vw[i].w > 0.0) { + tmp[nkeep] = vw[i]; + totw_keep += vw[i].w; + nkeep++; + } + } + if (nkeep == 0 || totw_keep <= 0.0) { + break; + } + + qsort(tmp, (size_t)nkeep, sizeof(valweight), cmp_valweight); + med = weighted_median_sorted(tmp, nkeep, totw_keep); + for (i = 0; i < nkeep; i++) { + tmp[i].v = fabs(tmp[i].v - med); + } + qsort(tmp, (size_t)nkeep, sizeof(valweight), cmp_valweight); + madv = weighted_median_sorted(tmp, nkeep, totw_keep); + sig = madv * mad_scale; + if (sig <= 0.0) { + break; + } + lo = med - clip_sigma * sig; + hi = med + clip_sigma * sig; + changed = 0; + for (i = 0; i < nvals; i++) { + if (keep[i] && (vw[i].v < lo || vw[i].v > hi)) { + keep[i] = 0; + changed = 1; + } + } + if (!changed) { + break; + } + } + + wsum = 0.0; + vsum = 0.0; + for (i = 0; i < nvals; i++) { + if (keep[i] && vw[i].w > 0.0) { + wsum += vw[i].w; + vsum += vw[i].w * vw[i].v; + } + } + *mean_clip = (wsum > 0.0) ? (vsum / wsum) : NAN; + + /* weighted median */ + qsort(vw, (size_t)nvals, sizeof(valweight), cmp_valweight); + *median = weighted_median_sorted(vw, nvals, totw); + + /* weighted MAD (unscaled) */ + for (i = 0; i < nvals; i++) { + vw[i].v = fabs(vw[i].v - *median); + } + qsort(vw, (size_t)nvals, sizeof(valweight), cmp_valweight); + madv = weighted_median_sorted(vw, nvals, totw); + *mad_std = madv * mad_scale; + +exit: + if (keep) { + free(keep); + } + if (tmp) { + free(tmp); + } + if (vw) { + free(vw); + } + + return status; +} + /*****************************************************************************/ /* elliptical annulus aperture */ diff --git a/src/sep.h b/src/sep.h index 746e8e5..5a8e102 100644 --- a/src/sep.h +++ b/src/sep.h @@ -306,6 +306,31 @@ SEP_API int sep_sum_circann( short * flag ); +/* sep_stats_circann() + * + * Compute statistics within a circular annulus. The MAD output is scaled + * by 1.4826, providing a robust estimate of the standard deviation for + * Gaussian data. + */ +SEP_API int sep_stats_circann( + const sep_image * image, + double x, + double y, + double rin, + double rout, + int id, + int subpix, + short inflags, + double clip_sigma, + int clip_iters, + double * mean, + double * std, + double * median, + double * mad_std, + double * mean_clip, + short * flag +); + SEP_API int sep_sum_ellipse( const sep_image * image, double x, diff --git a/test.py b/test.py index f507bab..4642374 100755 --- a/test.py +++ b/test.py @@ -784,6 +784,71 @@ def test_apertures_exact(): assert_allclose(flux, np.pi * ratio * (rout**2 - r**2)) +def _sigma_clip_mean(values, sigma=3.0, maxiters=5): + mask = np.ones(values.shape, dtype=bool) + for _ in range(maxiters): + vals = values[mask] + if vals.size == 0: + return np.nan + med = np.median(vals) + mad = np.median(np.abs(vals - med)) + std = 1.4826 * mad + if std == 0.0: + return vals.mean() + lo = med - sigma * std + hi = med + sigma * std + newmask = mask & (values >= lo) & (values <= hi) + if newmask.sum() == mask.sum(): + return vals.mean() + mask = newmask + vals = values[mask] + return vals.mean() if vals.size else np.nan + + +def test_stats_circann_flat(): + data = np.ones(data_shape) + rin, rout = 3.0, 6.0 + mean, std, med, mad_std, mean_clip, flag = sep.stats_circann( + data, x, y, rin, rout, subpix=1 + ) + assert_allclose(mean, 1.0) + assert_allclose(std, 0.0) + assert_allclose(med, 1.0) + assert_allclose(mad_std, 0.0) + assert_allclose(mean_clip, 1.0) + + +def test_stats_circann_matches_numpy_subpix1(): + rng = np.random.default_rng(12345) + data = rng.normal(size=(64, 64)) + x0 = np.array([30.4]) + y0 = np.array([31.2]) + rin, rout = 5.0, 8.0 + + mean, std, med, mad_std, mean_clip, flag = sep.stats_circann( + data, x0, y0, rin, rout, subpix=1 + ) + + yy, xx = np.indices(data.shape) + rpix2 = (xx - x0[0]) ** 2 + (yy - y0[0]) ** 2 + mask = (rpix2 >= rin**2) & (rpix2 < rout**2) + vals = data[mask] + + exp_mean = vals.mean() + exp_std = vals.std() + exp_med = np.median(vals) + exp_mad = np.median(np.abs(vals - exp_med)) + exp_mad_std = 1.4826 * exp_mad # robust std estimate for Gaussian data + exp_clip = _sigma_clip_mean(vals) + + assert_allclose(mean[0], exp_mean, rtol=1.0e-10, atol=1.0e-8) + assert_allclose(std[0], exp_std, rtol=1.0e-10, atol=1.0e-8) + assert_allclose(med[0], exp_med, rtol=1.0e-10, atol=1.0e-8) + assert_allclose(mad_std[0], exp_mad_std, rtol=1.0e-10, atol=1.0e-7) + assert_allclose(mean_clip[0], exp_clip, rtol=1.0e-10, atol=1.0e-8) + assert flag[0] == 0 + + def test_aperture_bkgann_overlapping(): """ Test bkgann functionality in circular & elliptical apertures. From 07462a3947abf5963e7f33438fd8854dbd8b0a82 Mon Sep 17 00:00:00 2001 From: karpov-sv Date: Mon, 2 Feb 2026 15:47:44 +0100 Subject: [PATCH 2/3] Couple of minor fixes --- sep.pyx | 6 +++--- src/aperture.i | 19 +++++++++++++++---- 2 files changed, 18 insertions(+), 7 deletions(-) diff --git a/sep.pyx b/sep.pyx index be241ec..b26dbcf 100644 --- a/sep.pyx +++ b/sep.pyx @@ -1092,7 +1092,7 @@ def sum_circann(np.ndarray data not None, x, y, rin, rout, var=None, err=None, gain=None, np.ndarray mask=None, double maskthresh=0.0, seg_id=None, np.ndarray segmap=None, int subpix=5): - """sum_circann(data, x, y, rin, rout, err=None, var=None, mask=None, + """sum_circann(data, x, y, rin, rout, var=None, err=None, mask=None, maskthresh=0.0, seg_id=None, segmap=None, gain=None, subpix=5) @@ -1111,8 +1111,8 @@ def sum_circann(np.ndarray data not None, x, y, rin, rout, element of the array. These inputs obey numpy broadcasting rules. It is required that ``rout >= rin >= 0.0``. - err, var : float or ndarray - Error *or* variance (specify at most one). + var, err : float or ndarray + Variance *or* error (specify at most one). mask : `~numpy.ndarray`, optional Mask array. If supplied, a given pixel is masked if its value diff --git a/src/aperture.i b/src/aperture.i index 53bd089..ac313b1 100644 --- a/src/aperture.i +++ b/src/aperture.i @@ -46,9 +46,15 @@ int APER_NAME( errort = im->noise; *flag = 0; varpix = 0.0; - scale = 1.0 / subpix; - scale2 = scale * scale; - offset = 0.5 * (scale - 1.0); + if (subpix > 0) { + scale = 1.0 / subpix; + scale2 = scale * scale; + offset = 0.5 * (scale - 1.0); + } else { + scale = 0.0; + scale2 = 0.0; + offset = 0.0; + } errisarray = 0; errisstd = 0; @@ -185,7 +191,12 @@ int APER_NAME( /* correct for masked values */ if (im->mask) { - if (inflag & SEP_MASK_IGNORE) { + if (totarea > 0.0 && maskarea >= totarea) { + *flag |= SEP_APER_ALLMASKED; + tv = 0.0; + sigtv = 0.0; + totarea = 0.0; + } else if (inflag & SEP_MASK_IGNORE) { totarea -= maskarea; } else { tv *= (tmp = totarea / (totarea - maskarea)); From c0523d3c9d8334164527de99e4a3b5c27c002a84 Mon Sep 17 00:00:00 2001 From: karpov-sv Date: Tue, 3 Feb 2026 11:23:47 +0100 Subject: [PATCH 3/3] sep_stats_circann: add area and sumerr outputs Compute annulus area and sumerr in one pass. Update the C API and Python wrapper. Document the new outputs. --- docs/apertures.rst | 4 +++ docs/changelogs/changes_to_c_api.rst | 5 ++++ docs/reference.rst | 1 + sep.pyx | 6 +++- src/aperture.c | 44 ++++++++++++++++++++++++++-- src/sep.h | 2 ++ 6 files changed, 58 insertions(+), 4 deletions(-) diff --git a/docs/apertures.rst b/docs/apertures.rst index 8fa7e30..fefcd27 100644 --- a/docs/apertures.rst +++ b/docs/apertures.rst @@ -12,6 +12,10 @@ Function Sums data within... `sep.sum_ellipann` elliptical annulus/annuli ================== ========================= + +The `~sep.stats_circann` function computes statistics in circular annuli, +including a sigma-clipped mean and robust scatter estimates. + The follow examples demonstrate options for circular aperture photometry. The other functions behave similarly. diff --git a/docs/changelogs/changes_to_c_api.rst b/docs/changelogs/changes_to_c_api.rst index 8e18e72..cbc457d 100644 --- a/docs/changelogs/changes_to_c_api.rst +++ b/docs/changelogs/changes_to_c_api.rst @@ -74,6 +74,11 @@ All changes here are transparent to users of the Python interface. - The type of ``n`` has changed from ``int`` to ``int64_t``. + +.. c:function:: int sep_stats_circann() + + - Added ``area`` and ``sumerr`` outputs. + .. c:function:: int sep_flux_radius() - The type of ``n`` has changed from ``int`` to ``int64_t``. diff --git a/docs/reference.rst b/docs/reference.rst index 07a114b..a0fb857 100644 --- a/docs/reference.rst +++ b/docs/reference.rst @@ -17,6 +17,7 @@ Reference/API sep.sum_circle sep.sum_circann + sep.stats_circann sep.sum_ellipse sep.sum_ellipann diff --git a/sep.pyx b/sep.pyx index b26dbcf..1379a05 100644 --- a/sep.pyx +++ b/sep.pyx @@ -170,7 +170,8 @@ cdef extern from "sep.h": int id, int subpix, short inflags, double clip_sigma, int clip_iters, double *mean, double *std, double *median, - double *mad_std, double *mean_clip, short *flag) + double *mad_std, double *mean_clip, + double *area, double *sumerr, short *flag) int sep_sum_ellipse(const sep_image *image, double x, double y, double a, double b, double theta, @@ -1311,6 +1312,7 @@ def stats_circann(np.ndarray data not None, x, y, rin, rout, """ cdef int status + cdef double area1, sumerr1 cdef np.broadcast it cdef sep_image im @@ -1366,6 +1368,8 @@ def stats_circann(np.ndarray data not None, x, y, rin, rout, np.PyArray_MultiIter_DATA(it, 7), np.PyArray_MultiIter_DATA(it, 8), np.PyArray_MultiIter_DATA(it, 9), + &area1, + &sumerr1, np.PyArray_MultiIter_DATA(it, 10)) _assert_ok(status) diff --git a/src/aperture.c b/src/aperture.c index b9883c6..f1ae85d 100644 --- a/src/aperture.c +++ b/src/aperture.c @@ -346,19 +346,23 @@ int sep_stats_circann( double * median, double * mad_std, double * mean_clip, + double * area, + double * sumerr, short * flag ) { PIXTYPE pix; double dx, dy, dx1, dy2, offset, scale, scale2, rpix2, overlap; double rin2, rin_in2, rin_out2, rout2, rout_in2, rout_out2; double totw, sumw, varw, diff, wsum, vsum; + double tv, sigtv, varpix, varpix_const; double med, madv, sig, totw_keep; double lo, hi; int64_t ix, iy, i, xmin, xmax, ymin, ymax, sx, sy, pos, size, msize, ssize, nkeep; - int64_t nvals, cap; + int64_t nvals, cap, esize; int status, ismasked, changed, iter; - const BYTE *datat, *maskt, *segt; - converter convert, mconvert, sconvert; + short errisarray, errisstd; + const BYTE *datat, *errort, *maskt, *segt; + converter convert, econvert = NULL, mconvert, sconvert; valweight *vw = NULL; valweight *tmp = NULL; char *keep = NULL; @@ -383,6 +387,8 @@ int sep_stats_circann( *median = NAN; *mad_std = NAN; *mean_clip = NAN; + *area = 0.0; + *sumerr = 0.0; rin2 = rin * rin; rout2 = rout * rout; @@ -403,6 +409,11 @@ int sep_stats_circann( segt = NULL; msize = 0; ssize = 0; + errort = im->noise; + esize = 0; + errisarray = 0; + errisstd = 0; + varpix_const = 0.0; /* get data converter(s) for input array(s) */ if ((status = get_converter(im->dtype, &convert, &size))) { @@ -429,12 +440,17 @@ int sep_stats_circann( nvals = 0; totw = 0.0; + tv = 0.0; + sigtv = 0.0; /* loop over rows in the box */ for (iy = ymin; iy < ymax; iy++) { /* set pointers to the start of this row */ pos = (iy % im->h) * im->w + xmin; datat = MSVC_VOID_CAST im->data + pos * size; + if (errisarray) { + errort = MSVC_VOID_CAST im->noise + pos * esize; + } if (im->mask) { maskt = MSVC_VOID_CAST im->mask + pos * msize; } @@ -506,6 +522,16 @@ int sep_stats_circann( vw[nvals].v = pix; vw[nvals].w = overlap; totw += overlap; + tv += pix * overlap; + if (errisarray) { + varpix = econvert(errort); + if (errisstd) { + varpix *= varpix; + } + } else { + varpix = varpix_const; + } + sigtv += varpix * overlap; nvals++; } } @@ -513,6 +539,9 @@ int sep_stats_circann( /* increment pointers by one element */ datat += size; + if (errisarray) { + errort += esize; + } maskt += msize; segt += ssize; } @@ -611,6 +640,15 @@ int sep_stats_circann( madv = weighted_median_sorted(vw, nvals, totw); *mad_std = madv * mad_scale; + if (im->gain > 0.0 && tv > 0.0) { + sigtv += tv / im->gain; + } + if (sigtv < 0.0) { + sigtv = 0.0; + } + *sumerr = sqrt(sigtv); + *area = totw; + exit: if (keep) { free(keep); diff --git a/src/sep.h b/src/sep.h index 5a8e102..131c2fb 100644 --- a/src/sep.h +++ b/src/sep.h @@ -328,6 +328,8 @@ SEP_API int sep_stats_circann( double * median, double * mad_std, double * mean_clip, + double * area, + double * sumerr, short * flag );