import json, numpy as np
from PIL import Image, ImageDraw
R = json.load(open("rhino.json"))
img = np.asarray(Image.open("base100-1.png").convert("RGB")).astype(int)
H, W, _ = img.shape
r, g, b = img[..., 0], img[..., 1], img[..., 2]
bleu = ((b - r) > 30) & (b > 150)                        # eau
vert = ((g - r) > 6) & ((g - b) > 25) & (g > 170)       # limites de parcelles (vert pâle)
cible = (bleu * 2.0 + vert * 1.0).astype(np.float32)
print("pixels bleus", bleu.sum(), "verts", vert.sum())
eau = [c for c in R["courbes"] if c["l"] == "B_EAU noues" and c["f"]]
parc = [c for c in R["courbes"] if c["l"] == "PARCELLES"]
allp = np.array([p for c in parc for p in c["p"]])
x0, y1 = allp[:, 0].min(), allp[:, 1].max()
def raster(s, shape):
    im = Image.new("F", (shape[1], shape[0]), 0); d = ImageDraw.Draw(im)
    for c in eau:
        d.polygon([((x - x0) * s, (y1 - y) * s) for x, y in c["p"]], fill=2.0)
    for c in parc:
        d.line([((x - x0) * s, (y1 - y) * s) for x, y in c["p"]], fill=1.0, width=max(1, int(s * 1.2)))
    return np.asarray(im, dtype=np.float32)
def correl(a, bimg):
    P = (2 * H, 2 * W)
    Fa = np.fft.rfft2(a, P); Fb = np.fft.rfft2(bimg, P)
    c = np.fft.irfft2(Fa * np.conj(Fb), P)
    k = np.unravel_index(np.argmax(c), c.shape)
    return c[k], k
best = None
for s in np.arange(2.0, 9.0, 0.05):
    rs = raster(s, (H, W))
    v, k = correl(cible, rs)
    v = v / (np.sqrt((rs ** 2).sum()) + 1e-6)
    if best is None or v > best[0]: best = (v, s, k)
v, s, k = best
dy = k[0] if k[0] < H else k[0] - 2 * H
dx = k[1] if k[1] < W else k[1] - 2 * W
print("grossier : s=%.3f px/m  dx=%d dy=%d  score=%.1f" % (s, dx, dy, v))
# affinage
best2 = None
for s2 in np.arange(s - 0.06, s + 0.06, 0.005):
    rs = raster(s2, (H, W)); v2, k2 = correl(cible, rs); v2 /= (np.sqrt((rs ** 2).sum()) + 1e-6)
    if best2 is None or v2 > best2[0]: best2 = (v2, s2, k2)
v, s, k = best2
dy = k[0] if k[0] < H else k[0] - 2 * H
dx = k[1] if k[1] < W else k[1] - 2 * W
print("fin : s=%.4f px/m à 100 dpi  dx=%d dy=%d" % (s, dx, dy))
# transformation mètres -> pixels (100 dpi) : px = (x - x0)*s + dx ; py = (y1 - y)*s + dy
T = {"dpi": 100, "s": float(s), "x0": float(x0), "y1": float(y1), "dx": int(dx), "dy": int(dy)}
json.dump(T, open("transfo.json", "w"))
# contrôle visuel : géométrie Rhino en rouge par-dessus le dessin
ov = Image.open("base100-1.png").convert("RGB"); d = ImageDraw.Draw(ov)
f = lambda x, y: ((x - x0) * s + dx, (y1 - y) * s + dy)
for c in parc: d.line([f(*p) for p in c["p"]], fill=(220, 30, 30), width=1)
for c in eau: d.line([f(*p) for p in c["p"]] + [f(*c["p"][0])], fill=(200, 0, 200), width=2)
ov.save("controle-recalage.png")