"""Tvorba inverzní (negativní) jednodílné formy z korigovaného reliéfu.

Forma je těleso  { (x,y,z) : h(x,y) <= z <= O(x,y) }  nad doménou lemu, kde
h je korigovaný reliéf (dutina formy) a O je jeho **skutečný ekvidistantní
offset** o tloušťku stěny t.

Offset se počítá jako morfologická dilatace tělesa koulí o poloměru t,
která je pro výškové pole opět výškovým polem:

    O(p) = max_{|q-p| <= t} ( h(q) + sqrt(t^2 - |q-p|^2) )

Tím je zaručena tloušťka stěny >= t měřená kolmo k povrchu (nikoli jen
svisle) — svislé stěny reliéfu tak nedostanou tenkou blánu.
"""
from __future__ import annotations

import numpy as np

from .draft import cone_dilate
from .raster import Grid


def _ball_structure(res: float, t: float, conservative: bool = True):
    """Diskrétní strukturní funkce koule sqrt(t^2 - d^2).

    Při conservative=True se vzdálenost zmenší o poloviční úhlopříčku pixelu,
    takže diskrétní jádro leží *nad* spojitým — výsledný plášť je nanejvýš
    nepatrně silnější, nikdy tenčí.
    """
    r = int(np.ceil(t / res)) + (1 if conservative else 0)
    ax = np.arange(-r, r + 1) * res
    dx, dy = np.meshgrid(ax, ax, indexing="xy")
    d = np.hypot(dx, dy)
    if conservative:
        d = np.maximum(d - 0.70711 * res, 0.0)
    fp = d <= t + 1e-12
    st = np.zeros_like(d)
    st[fp] = np.sqrt(np.maximum(t * t - d[fp] ** 2, 0.0))
    return fp, st, r


def _grey_ball(h: np.ndarray, res: float, t: float,
               conservative: bool = True) -> np.ndarray:
    from scipy.ndimage import grey_dilation
    fp, st, _ = _ball_structure(res, t, conservative)
    return grey_dilation(h, footprint=fp, structure=st, mode="nearest")


def _maxpool(h: np.ndarray, c: int) -> np.ndarray:
    ny, nx = h.shape
    py = (-ny) % c
    px = (-nx) % c
    if py or px:
        h = np.pad(h, ((0, py), (0, px)), mode="edge")
    ny2, nx2 = h.shape
    return h.reshape(ny2 // c, c, nx2 // c, c).max(axis=(1, 3))


def _upsample(coarse: np.ndarray, shape, c: int) -> np.ndarray:
    from scipy.ndimage import map_coordinates
    ny, nx = shape
    jj = (np.arange(ny) + 0.5) / c - 0.5
    ii = (np.arange(nx) + 0.5) / c - 0.5
    J, I = np.meshgrid(jj, ii, indexing="ij")
    return map_coordinates(coarse, [J, I], order=1, mode="nearest")


def ball_offset(h: np.ndarray, res: float, t: float,
                method: str = "auto", target_r_px: int = 16) -> np.ndarray:
    """Ekvidistantní offset výškového pole o t = dilatace koulí o poloměru t.

        O(p) = max_{|q-p| <= t} ( h(q) + sqrt(t^2 - |q-p|^2) )

    method='exact'  — přesná morfologická dilatace v plném rozlišení
    method='fast'   — dilatace na zhrubené mřížce (h se zhrubuje **maximem**,
                      výsledek je proto vždy >= přesný offset, tj. plášť je
                      nanejvýš o zlomek mm silnější, nikdy tenčí) s následnou
                      bilineární interpolací zpět.
    method='auto'   — 'exact' pro malé poloměry v pixelech, jinak 'fast'.
    """
    if t <= 0:
        return h.copy()
    r = int(np.ceil(t / res))
    if method == "auto":
        method = "exact" if r <= target_r_px else "fast"

    try:
        if method == "exact" or r <= target_r_px:
            return np.maximum(_grey_ball(h, res, t), h + t)
        c = int(np.ceil(r / float(target_r_px)))
        # kompenzace zhrubení: vzorkovací posun + bilineární interpolace
        t_c = t + 0.5 * c * res
        coarse = _maxpool(h, c)
        off = _grey_ball(coarse, res * c, t_c)
        up = _upsample(off, h.shape, c)
        return np.maximum(up, h + t)
    except Exception:
        return h + t


def footprint_mask(h: np.ndarray, eps: float = 1e-9) -> np.ndarray:
    from .meshing import fill_holes
    return fill_holes(h > eps)


def rim_domain(fp: np.ndarray, res: float, rim_mm: float) -> np.ndarray:
    """Doména formy = půdorys reliéfu roztažený o šířku lemu."""
    if rim_mm <= 0:
        return fp.copy()
    try:
        from scipy.ndimage import distance_transform_edt
        dist = distance_transform_edt(~fp, sampling=(res, res))
        return fp | (dist <= rim_mm)
    except Exception:
        # nouzová varianta: kuželová dilatace indikátoru
        ind = np.where(fp, 0.0, -1e9)
        grown = cone_dilate(ind, res, 1.0 / max(rim_mm, 1e-9), rounds=2)
        return grown > -1.0


def smooth_outline(dom: np.ndarray, res: float, radius_mm: float) -> np.ndarray:
    """Morfologické uzavření obrysu — zahladí úzké vruby a ostré vnitřní rohy.

    Uzavření (dilatace + eroze týmž poloměrem) může doménu jen zvětšit, nikdy
    zmenšit, takže lem nikde neztenčí pod zadanou šířku.
    """
    if radius_mm <= 0:
        return dom
    try:
        from scipy.ndimage import distance_transform_edt
        grown = distance_transform_edt(~dom, sampling=(res, res)) <= radius_mm
        shrunk = distance_transform_edt(grown, sampling=(res, res)) >= radius_mm
        return dom | shrunk
    except Exception:
        return dom


def build_mold(grid: Grid, h: np.ndarray, wall_mm: float = 3.0,
               rim_mm: float = 5.0, back: str = "shell",
               offset_method: str = "auto", rim_smooth_mm: float = -1.0):
    """Vrátí (z_bottom, z_top, cell_mask, info) pro meshování formy."""
    from .meshing import cells_from_nodes_all, fill_holes

    fp = footprint_mask(h)
    dom = fill_holes(rim_domain(fp, grid.res, rim_mm))
    r_smooth = rim_mm if rim_smooth_mm is not None and rim_smooth_mm < 0 \
        else float(rim_smooth_mm)
    dom = fill_holes(smooth_outline(dom, grid.res, r_smooth))

    if back == "flat":
        top = np.full_like(h, float(h.max()) + wall_mm)
    else:
        top = ball_offset(h, grid.res, wall_mm, offset_method)
        top = np.maximum(top, h + wall_mm)

    bottom = np.maximum(h, 0.0)
    cm = cells_from_nodes_all(dom)
    if not cm.any():
        cm = cells_from_nodes_all(fp)

    info = {
        "footprint_area_mm2": float(fp.sum()) * grid.cell_area,
        "domain_area_mm2": float(dom.sum()) * grid.cell_area,
        "mold_height_mm": float(top[dom].max()) if dom.any() else 0.0,
        "min_wall_mm": float(np.min((top - bottom)[dom])) if dom.any() else 0.0,
    }
    return bottom, top, cm, info
