import os
import glob
import numpy as np
from PIL import Image

def _ensure_float(array):
    """Return a contiguous float64 copy of the provided array."""
    return np.asarray(array, dtype=np.float64, order="C")


def read_image(folder, name, scale=1, gray=False):
    for filename in glob.glob(folder + "/*"):
        if os.path.splitext(os.path.basename(filename))[0] == name:
            break
    img = Image.open(os.path.join(filename))
    if scale != 1:
        w, h = img.size
        img = img.resize((int(w * scale), int(h * scale)), Image.BICUBIC)
    if gray:
        img = img.convert("L")
    img = np.array(img)
    if len(img.shape) == 3:
        img = img[..., :3]
    return img.astype(np.float64) / 255 # only first 3

def rgb2gray(img):
    r, g, b = img[..., 0], img[..., 1], img[..., 2]
    return 0.2126 * r + 0.7152 * g + 0.0722 * b

def rgb2lum(img):
    """
    Simple luminance helper used by Perona–Malik code.
    If img is RGB, return standard luma; if already single-channel, return as-is.
    """
    img = np.asarray(img)
    if img.ndim == 3 and img.shape[-1] == 3:
        return rgb2gray(img)
    return img

def forward_gradients(img):
    """
    Forward finite-difference gradients with zero-Neumann boundary conditions.
    """
    img = _ensure_float(img)
    dx = np.zeros_like(img)
    dy = np.zeros_like(img)
    dx[:, :-1] = img[:, 1:] - img[:, :-1]
    dy[:-1, :] = img[1:, :] - img[:-1, :]
    return dx, dy


def gradient_magnitude(img):
    """Compute sqrt(dx^2 + dy^2) using forward differences."""
    dx, dy = forward_gradients(img)
    return np.sqrt(dx * dx + dy * dy + 1e-12)


def perona_malik_conductivity_from_guidance(guidance, kappa):
    """
    Build scalar Perona–Malik conductivities c_p = exp(-(|∇h| / κ)^2).
    """
    grad_mag = gradient_magnitude(guidance)
    scaled = grad_mag / max(kappa, 1e-8)
    return np.exp(-(scaled ** 2))


def edge_weights_from_conductivity(conductivity):
    """
    Convert per-pixel conductivities into horizontal/vertical edge weights.
    """
    c = _ensure_float(conductivity)
    wx = 0.5 * (c[:, :-1] + c[:, 1:])
    wy = 0.5 * (c[:-1, :] + c[1:, :])
    return wx, wy


def build_mask_index(mask):
    """
    Return an index map such that idx[y, x] = position inside Ω_M or -1 outside.
    """
    idx = -np.ones_like(mask, dtype=np.int32)
    idx[mask.astype(bool)] = np.arange(mask.sum(), dtype=np.int32)
    return idx

