Documentation
¶
Overview ¶
Package stats provides generic, NaN-aware statistics functions for astronomical image data. Every function operates on slices of any Numeric type and uses float64 accumulators internally for precision.
The core algorithms are extracted from the library's existing cross-validated implementations:
- Percentile/Median/MAD: ported from Siril/RawTherapee's findMinMaxPercentile (histogram-based interpolated percentile)
- MeanStdev: ported from cfitsio's FnMeanSigma
NaN handling: for float32/float64 inputs, NaN values are silently skipped in all computations. Integer types cannot be NaN, so the generic code fast-paths them with no NaN checks.
Index ¶
- Constants
- func FilterNonZero[T Numeric](data []T) []T
- func MAD[T Numeric](data []T) float64
- func MADWithMedian[T Numeric](data []T) (median, mad float64)
- func MADWithMedianBuf[T Numeric](data []T, histo []uint32, absDev []float64) (median, mad float64)
- func Mean[T Numeric](data []T) float64
- func MeanStdev[T Numeric](data []T) (mean, stdev float64)
- func MeanStdevSiril(data []float32) (mean, stdev float32)
- func Median[T Numeric](data []T) float64
- func MedianBuf[T Numeric](data []T, histo []uint32) float64
- func MinMax[T Numeric](data []T) (min, max T)
- func Percentile[T Numeric](data []T, p float64) float64
- func PercentileBuf[T Numeric](data []T, p float64, histo []uint32) float64
- func QuickSelect[T Numeric](data []T, k int) T
- type CenterFunc
- type Histogram
- type Numeric
- type SigmaClipResult
Constants ¶
const MaxHistoSize = 65536
MaxHistoSize is the upper bound on histogram bins used by the percentile/median family. Callers using the *Buf zero-alloc variants should pre-allocate a histogram slice of this length (or len(data), whichever is smaller).
Variables ¶
This section is empty.
Functions ¶
func FilterNonZero ¶
func FilterNonZero[T Numeric](data []T) []T
FilterNonZero returns a new slice containing only the non-zero, non-NaN elements of data. This matches the astronomical convention for "blank pixel" exclusion — Siril's reassign_to_non_null_data_float uses the same logic before computing median and MAD.
For integer types, excludes exact zeros (NaN is not possible). For float types, excludes both exact 0.0 and NaN.
If all elements are zero/NaN, returns the original slice as a fallback (avoids returning empty data to downstream statistics).
func MAD ¶
MAD returns the Median Absolute Deviation: median(|xi - median(x)|). This is the raw MAD — caller multiplies by 1.4826 for Gaussian- equivalent sigma if needed.
Allocates two slices on every call (histogram + absolute deviations). For hot paths, use MADWithMedianBuf with caller-supplied scratch.
func MADWithMedian ¶
MADWithMedian returns the median and the MAD, avoiding double- computing the median when the caller needs both.
Allocates two slices on every call. For hot paths, use MADWithMedianBuf.
func MADWithMedianBuf ¶ added in v1.1.0
MADWithMedianBuf is the zero-alloc variant of MADWithMedian. Caller supplies two scratch slices:
- histo: must have len >= min(len(data), MaxHistoSize); used for both the median pass and the MAD-of-deviations pass. Zeroed and overwritten internally.
- absDev: must have cap >= len(data); used as scratch for the |xi - median| array. Length is reset to 0 on entry, then grown by append up to the non-NaN count of data.
Panics if either buffer is too small.
Note: absDev is []float64 (not generic []T) so a single buffer works for any T. The deviations are computed as float64 regardless of input type, matching the existing MADWithMedian.
func Mean ¶
Mean returns the arithmetic mean of data, skipping NaN. For empty slices (or all-NaN), returns 0.
func MeanStdev ¶
MeanStdev returns the arithmetic mean and sample standard deviation (Bessel-corrected, N-1 denominator) of data, skipping NaN. Uses float64 accumulators. For empty or single-element slices, stdev is 0.
The N-1 denominator matches Siril's siril_stats_float_sd, which is the oracle we cross-validate against.
func MeanStdevSiril ¶ added in v1.1.0
MeanStdevSiril computes the mean and sample standard deviation (N-1) of data using the two-pass algorithm from Siril's siril_stats_float_sd (statistics_float.c). Returns float32 to preserve rejection-boundary precision — silently promoting to float64 produces ULP-level differences at sigma-clip boundaries that affect pixel rejection decisions on real astronomical data.
Algorithm (bit-for-bit port of the C source):
- First pass: accumulate sum in float64.
- Truncate the mean to float32: mean = float32(sum / n).
- Second pass: subtract in float32 (data[i] - mean), square in float32, accumulate in float64.
- Return float32(sqrt(acc / (n-1))).
The float32 mean truncation and float32 deviation arithmetic are the load-bearing parts: the same loops in float64 produce different ULP- level results that change pixel rejection decisions. The fits- processing rejection pipeline traced a visible artifact (red hotspot in stacked output) to this precision gap and validated the two-pass- with-truncation algorithm against Siril's actual rejection output on real Rosette Nebula data.
NaN handling: NONE — matches Siril. NaN inputs propagate by IEEE rules and the result is NaN in both outputs. Callers must pre-filter NaN if needed (cf. Siril's reassign_to_non_null_data_float).
Edge cases (also Siril behavior):
- n=0: returns (NaN, -0). Mean is 0/0 = NaN; the deviation accumulator is 0 and the (n-1)=-1 division yields -0.
- n=1: returns (data[0], NaN). Mean is data[0]; the single deviation is 0 and the (n-1)=0 division yields NaN.
If you want a friendly default (e.g. (0, 0) for empty input), use the generic MeanStdev — but be aware its outputs differ from this function at the ULP level even for well-formed inputs.
Source: Siril src/algos/statistics_float.c siril_stats_float_sd.
func Median ¶
Median returns the median of data (the 50th percentile).
Allocates a histogram on every call. For hot paths, use MedianBuf with a caller-supplied scratch buffer.
func MedianBuf ¶ added in v1.1.0
MedianBuf is the zero-alloc variant of Median. The histo slice must have len >= min(len(data), MaxHistoSize); panics otherwise. histo is zeroed and overwritten on each call.
func MinMax ¶
func MinMax[T Numeric](data []T) (min, max T)
MinMax returns the minimum and maximum of data, skipping NaN. For empty slices, returns (0, 0).
func Percentile ¶
Percentile returns the p-th percentile (p in [0, 1], clamped) of data using a histogram-based interpolated algorithm. This is a line-for-line port of RawTherapee/Siril's findMinMaxPercentile (rt_algo.cc:167), made generic over any Numeric type.
The algorithm:
- Find min/max of the non-NaN data
- Build a histogram with min(65536, n) bins spanning [min, max]
- Walk the CDF to the percentile
- Interpolate between the two bins straddling the percentile
- Convert back to the original value range
For empty slices (or all-NaN), returns 0.
Allocates a histogram slice on every call. For hot paths, use PercentileBuf with a caller-supplied scratch buffer.
func PercentileBuf ¶ added in v1.1.0
PercentileBuf is the zero-alloc variant of Percentile. The histo slice must have len >= min(len(data), MaxHistoSize); panics otherwise. histo is zeroed and overwritten on each call, so callers can reuse the same buffer across many calls without resetting it.
Typical usage: pre-allocate one MaxHistoSize histogram per goroutine and reuse it for every Percentile/Median/MAD call:
var histo = make([]uint32, stats.MaxHistoSize)
for _, frame := range frames {
med := stats.PercentileBuf(frame, 0.5, histo)
// ...
}
func QuickSelect ¶ added in v1.1.0
QuickSelect rearranges data so that data[k] contains the k-th smallest element (0-indexed). Returns data[k]. Uses the Hoare-style two-pointer partition with the middle element as pivot. O(n) average, O(n²) worst case.
This is a generic line-for-line port of Siril's quickmedian_float (sorting.c) — same partition, same pivot strategy, same loop structure. The fits-processing rejection pipeline already has a float32-specific copy (stack/reject.go quickselectF32); this is the library version, generic over any Numeric type. Behavioral equivalence with the float32 version is intentional, so callers can drop the local copy and use this without changing rejection decisions.
Unlike Median, this computes the exact order statistic — no histogram quantization or sub-bin interpolation. For small n (under ~50 elements) this is both faster than Median and more accurate, since the histogram method's bin width is comparable to the data spacing in that regime. For large n, prefer MedianBuf (zero-alloc, O(n)) which is faster despite the histogram overhead.
IMPORTANT differences from Median for the same input:
- QuickSelect returns an exact element of data, not an interpolated value between two elements
- For even n, "the median" is conventionally the mean of the two middle elements; QuickSelect with k=n/2 returns just the upper of those two
- Therefore QuickSelect(data, len(data)/2) is NOT equivalent to Median(data); callers who need a drop-in replacement should compute (QuickSelect(data, n/2-1) + QuickSelect(data, n/2)) / 2 for even n
NaN handling: NONE. NaN comparisons always return false, which breaks partitioning — NaN values can end up anywhere in the output and corrupt the answer. Callers MUST pre-filter NaN before calling.
Mutation: data is rearranged in place. If callers need to preserve the original order, copy the slice first.
Panics if data is empty, or if k < 0 or k >= len(data).
Source: Siril src/algos/sorting.c quickmedian_float, mirrored at fits-processing/stack/reject.go quickselectF32.
Types ¶
type CenterFunc ¶
type CenterFunc int
CenterFunc selects the center estimator for sigma clipping.
const ( CenterMean CenterFunc = iota // Use arithmetic mean as center CenterMedian // Use median as center (more robust) )
type Histogram ¶
type Histogram struct {
Counts []int64 // bin counts, length = NBins
Edges []float64 // bin edges, length = NBins + 1
NBins int
Total int64 // sum of all counts
}
Histogram holds the result of a histogram computation.
func BuildHistogram ¶
BuildHistogram bins data into nbins equal-width bins spanning [min, max] of the non-NaN data. Returns a Histogram struct.
func (*Histogram) CDF ¶
CDF returns the cumulative distribution function as a float64 slice of length NBins. cdf[i] is the fraction of data at or below the upper edge of bin i.
func (*Histogram) Percentile ¶
Percentile returns the p-th percentile (p in [0, 1]) by interpolating the CDF of this histogram.
type Numeric ¶
type Numeric interface {
~uint8 | ~int8 | ~int16 | ~uint16 | ~int32 | ~uint32 | ~int64 | ~uint64 | ~float32 | ~float64
}
Numeric is the type constraint for all stats functions, matching the fits package's pixel type set.
type SigmaClipResult ¶
type SigmaClipResult struct {
Mean float64
Stdev float64
Median float64
NGood int // count of values that survived clipping
}
SigmaClipResult holds the output of iterative sigma clipping.
func SigmaClip ¶
func SigmaClip[T Numeric](data []T, kLow, kHigh float64, maxIter int, center CenterFunc) SigmaClipResult
SigmaClip performs iterative sigma clipping on data.
At each iteration, values outside [center - kLow*sigma, center + kHigh*sigma] are rejected. The center is recomputed from the surviving values using the chosen CenterFunc. Iteration stops when no more values are rejected or maxIter is reached (0 = iterate until stable, capped at 100).
Data is copied internally — the caller's slice is never modified. NaN values are excluded before clipping begins.
Source Files