From 69f61187e9f31e6be0eef31442e42d421c880209 Mon Sep 17 00:00:00 2001 From: Bojan Nikolic Date: Wed, 16 Jul 2025 15:52:33 +0000 Subject: [PATCH] Improve performance of iterative sigma clipped rms Sort the array first, then reuse for both median calculation and finding the clipped subset --- bdsf/functions.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/bdsf/functions.py b/bdsf/functions.py index a578fd4f..719ffa16 100644 --- a/bdsf/functions.py +++ b/bdsf/functions.py @@ -2135,6 +2135,7 @@ def bstat(indata, mask, kappa_npixbeam): unmasked = numpy.where(~msk_flat) skpix = skpix[unmasked] + skpix.sort() ct = skpix.size iter = 0 c1 = 1.0 @@ -2154,19 +2155,21 @@ def bstat(indata, mask, kappa_npixbeam): if kappa < 3.0: kappa = 3.0 lastct = ct - medval = numpy.median(skpix) + medval = skpix[npix//2] sig = numpy.std(skpix) - wsm = numpy.where(abs(skpix-medval) < kappa*sig) - ct = len(wsm[0]) + + wsm1, wsm2 = numpy.searchsorted(skpix, [medval - kappa*sig, + medval + kappa*sig]) + ct = wsm2 - wsm1 if ct > 0: - skpix = skpix[wsm] + skpix = skpix[wsm1:wsm2] c1 = abs(ct - lastct) c2 = converge_num * lastct iter += 1 mean = numpy.mean(skpix) - median = numpy.median(skpix) + median = skpix[skpix.size//2] sigma = numpy.std(skpix, ddof=1) mode = 2.5*median - 1.5*mean if sigma > 0.0: