binapprox.py

from astropy.io import fits
import numpy as np
import time
from statistics import median

# returns mean for each pixel position in files
def median_fits(files):
    start = time.perf_counter()
    data = np.array([ fits.open(f)[0].data for f in files ])
    _, h, w = data.shape
    r = np.zeros((h, w))
    for y in range(h):
        for x in range(w):
            r[y,x] = median(data[:,y,x])
    return r, time.perf_counter() - start, data.size * r.itemsize / 1024


def median_bins(values, B):
    mean    = np.mean(values)
    stdev   = np.std(values)
    min_val = mean - stdev
    max_val = mean + stdev

    left_bin = 0
    bins = np.zeros(B)
    bin_width = 2 * stdev / B

    for v in values:
        if v < min_val:
            left_bin += 1
        elif v < max_val:
            cur = int((v - min_val) / bin_width)
            bins[cur] += 1

    return mean, stdev, left_bin, bins

# Sum these counts until total >= (N + 1)/2. Remember to start from the ignore bin;
# Return the midpoint of the bin that exceeded (N + 1)/2.
def median_approx(values, B):
    mean, stdev, left_bin, bins = median_bins(values, B)

    n = len(values)
    mid = (n + 1) / 2

    count = left_bin
    for b, bin_count in enumerate(bins):
        count += bin_count
        if count >= mid:
            break

    bin_width = 2 * stdev / B
    median = mean - stdev + bin_width * (b + 0.5)

    return median

if __name__ == '__main__':
    print(median_bins([1, 1, 3, 2, 2, 6], 3))
    print(median_approx([1, 1, 3, 2, 2, 6], 3))

    print(median_bins([1, 5, 7, 7, 3, 6, 1, 1], 4))
    print(median_approx([1, 5, 7, 7, 3, 6, 1, 1], 4))