# Grouper: concentric circle radii by radial histograms of line pixels (plain Python). Prints centre refits and radii.
import numpy as np
from PIL import Image
import os, sys
HERE = os.path.dirname(os.path.abspath(__file__))
IM = np.asarray(Image.open(os.path.join(HERE, "..", "reference", "user_blueprint.jpg")).convert("RGB")).astype(float)
R = IM[..., 0]

def fit(cx, cy, r0, r1, thr=110):
    """least-squares circle through line pixels in the annulus r0..r1 (iterated twice)"""
    for _ in range(3):
        ys, xs = np.mgrid[int(cy - r1 - 2):int(cy + r1 + 3), int(cx - r1 - 2):int(cx + r1 + 3)]
        rr = np.hypot(xs - cx, ys - cy)
        m = (rr >= r0) & (rr <= r1) & (R[ys, xs] > thr)
        x, y = xs[m].astype(float), ys[m].astype(float)
        A = np.c_[2 * x, 2 * y, np.ones_like(x)]
        b = x * x + y * y
        (a, c, d), *_ = np.linalg.lstsq(A, b, rcond=None)
        cx, cy = a, c
        r = np.sqrt(d + a * a + c * c)
    return cx, cy, r, int(m.sum())

def radial(cx, cy, rmax, thr=110, step=0.5, sectors=None):
    ys, xs = np.mgrid[int(cy - rmax - 1):int(cy + rmax + 2), int(cx - rmax - 1):int(cx + rmax + 2)]
    rr = np.hypot(xs - cx, ys - cy)
    ang = np.degrees(np.arctan2(-(ys - cy), xs - cx)) % 360
    m = (rr <= rmax)
    if sectors:
        ms = np.zeros_like(m)
        for a0, a1 in sectors:
            ms |= (ang >= a0) & (ang <= a1)
        m &= ms
    w = np.clip(R[ys, xs] - 60, 0, None) * m
    bins = np.arange(0, rmax + step, step)
    h, _ = np.histogram(rr[m], bins=bins, weights=w[m])
    n, _ = np.histogram(rr[m], bins=bins)
    prof = h / np.maximum(n, 1)
    peaks = [(round(float(bins[i] + step / 2), 1), int(prof[i])) for i in range(1, len(prof) - 1)
             if prof[i] > prof[i - 1] and prof[i] >= prof[i + 1] and prof[i] > 45]
    return peaks

if __name__ == "__main__":
    a = list(map(float, sys.argv[1:]))
    cx, cy, rmax = a[:3]
    print("peaks", radial(cx, cy, rmax))