import base64
import math
from io import BytesIO

import cv2
import numpy as np
from numpy.lib.stride_tricks import sliding_window_view
from PIL import Image


LAPLACIAN_MASK = np.array([[0, -1, 0], [-1, 4, -1], [0, -1, 0]], dtype=np.float32)
SOBEL_GX = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]], dtype=np.float32)
SOBEL_GY = SOBEL_GX.T
ROBERTS_GX = np.array([[1, 0], [0, -1]], dtype=np.float32)
ROBERTS_GY = np.array([[0, 1], [-1, 0]], dtype=np.float32)


class ProcessingError(ValueError):
    pass


def ensure_uint8(image):
    return np.clip(image, 0, 255).astype(np.uint8)


def normalize_to_uint8(image):
    arr = image.astype(np.float32)
    min_value = float(np.min(arr))
    max_value = float(np.max(arr))
    if math.isclose(min_value, max_value):
        return np.zeros(arr.shape, dtype=np.uint8)
    return np.round((arr - min_value) * 255.0 / (max_value - min_value)).astype(np.uint8)


def require_odd(value, name="size", minimum=3):
    try:
        value = int(value)
    except (TypeError, ValueError) as exc:
        raise ProcessingError(f"{name} must be an odd integer.") from exc
    if value < minimum or value % 2 == 0:
        raise ProcessingError(f"{name} must be an odd integer >= {minimum}.")
    return value


def require_finite_positive(value, name):
    try:
        value = float(value)
    except (TypeError, ValueError) as exc:
        raise ProcessingError(f"{name} must be a finite positive number.") from exc
    if not np.isfinite(value) or value <= 0:
        raise ProcessingError(f"{name} must be a finite positive number.")
    return value


def to_gray(image):
    if image.ndim == 2:
        return image
    return cv2.cvtColor(image, cv2.COLOR_RGB2GRAY)


def gray_to_rgb(gray):
    return cv2.cvtColor(gray, cv2.COLOR_GRAY2RGB)


def histogram(image):
    gray = to_gray(image)
    counts = np.bincount(gray.ravel(), minlength=256).astype(np.float64)
    probabilities = counts / max(gray.size, 1)
    return probabilities.round(8).tolist()


def histogram_payload(image):
    gray = to_gray(image)
    payload = {"intensity": histogram(gray)}
    if image.ndim == 3:
        payload["r"] = (np.bincount(image[:, :, 0].ravel(), minlength=256).astype(np.float64) / image[:, :, 0].size).round(8).tolist()
        payload["g"] = (np.bincount(image[:, :, 1].ravel(), minlength=256).astype(np.float64) / image[:, :, 1].size).round(8).tolist()
        payload["b"] = (np.bincount(image[:, :, 2].ravel(), minlength=256).astype(np.float64) / image[:, :, 2].size).round(8).tolist()
    return payload


def image_to_data_url(image):
    pil_image = Image.fromarray(ensure_uint8(image))
    buffer = BytesIO()
    pil_image.save(buffer, format="PNG")
    payload = base64.b64encode(buffer.getvalue()).decode("ascii")
    return f"data:image/png;base64,{payload}"


def data_url_to_bytes(value):
    if "," in value:
        value = value.split(",", 1)[1]
    return base64.b64decode(value)


def decode_image(uploaded_file=None, base64_image=None):
    if uploaded_file is None and not base64_image:
        raise ProcessingError("Provide an image file or base64 image payload.")
    if uploaded_file is not None:
        raw = uploaded_file.read()
    else:
        raw = data_url_to_bytes(base64_image)

    image = Image.open(BytesIO(raw))
    image = image.convert("RGB")
    return np.array(image, dtype=np.uint8)


def negative(image, params):
    return 255 - image


def logarithmic(image, params):
    c = require_finite_positive(params.get("c", 1.0 / math.log(2.0)), "c")
    normalized = image.astype(np.float32) / 255.0
    transformed = c * np.log1p(normalized)
    return ensure_uint8(np.round(np.clip(transformed, 0.0, 1.0) * 255.0))


def gamma(image, params):
    gamma_value = require_finite_positive(params.get("gamma", 1.0), "gamma")
    c = require_finite_positive(params.get("c", 1.0), "c")
    normalized = image.astype(np.float32) / 255.0
    transformed = c * np.power(normalized, gamma_value)
    return ensure_uint8(np.round(np.clip(transformed, 0.0, 1.0) * 255.0))


def contrast_stretch(image, params):
    low = int(params.get("low", 0))
    high = int(params.get("high", 255))
    if low < 0 or high > 255 or low >= high:
        raise ProcessingError("Contrast stretch requires 0 <= low < high <= 255.")
    stretched = (image.astype(np.float32) - low) * (255.0 / (high - low))
    return ensure_uint8(np.round(stretched))


def gray_slice(image, params):
    start = int(params.get("start", 96))
    end = int(params.get("end", 160))
    if start < 0 or end > 255 or start > end:
        raise ProcessingError("Gray-level slicing requires 0 <= start <= end <= 255.")
    preserve = bool(params.get("preserve_background", True))
    highlight = np.array(params.get("highlight", [255, 64, 64]), dtype=np.uint8)
    if highlight.shape != (3,):
        raise ProcessingError("highlight must be an RGB triplet.")
    gray = to_gray(image)
    mask = (gray >= start) & (gray <= end)
    base = image.copy() if image.ndim == 3 else gray_to_rgb(gray if preserve else np.zeros_like(gray))
    if not preserve:
        base = np.zeros((*gray.shape, 3), dtype=np.uint8)
    base[mask] = highlight
    return base


def bit_plane(image, params):
    bit = int(params.get("bit", 7))
    if bit < 0 or bit > 7:
        raise ProcessingError("bit must be between 0 and 7.")
    plane = ((to_gray(image) >> bit) & 1) * 255
    return gray_to_rgb(plane.astype(np.uint8))


def histogram_equalization(image, params):
    gray = to_gray(image)
    counts = np.bincount(gray.ravel(), minlength=256)
    cdf = counts.cumsum().astype(np.float64)
    nonzero = cdf[cdf > 0]
    if nonzero.size == 0:
        return gray_to_rgb(gray)
    cdf_min = nonzero[0]
    denom = gray.size - cdf_min
    if denom <= 0:
        equalized = np.zeros_like(gray)
    else:
        lut = np.round((cdf - cdf_min) / denom * 255.0).clip(0, 255).astype(np.uint8)
        equalized = lut[gray]
    return gray_to_rgb(equalized)


def target_cdf_from_params(params):
    if "cdf" in params:
        cdf = np.array(params["cdf"], dtype=np.float64)
        if cdf.shape != (256,) or np.any(np.diff(cdf) < 0):
            raise ProcessingError("cdf must contain 256 non-decreasing values.")
        if cdf[-1] <= 0:
            raise ProcessingError("cdf must end with a positive value.")
        return cdf / cdf[-1]

    mode = params.get("target", "uniform")
    levels = np.arange(256, dtype=np.float64)
    if mode == "dark":
        pdf = np.exp(-levels / 64.0)
    elif mode == "bright":
        pdf = np.exp(-(255.0 - levels) / 64.0)
    elif mode == "bimodal":
        pdf = np.exp(-((levels - 72.0) ** 2) / (2 * 22.0**2)) + np.exp(-((levels - 190.0) ** 2) / (2 * 28.0**2))
    else:
        pdf = np.ones(256, dtype=np.float64)
    cdf = np.cumsum(pdf)
    return cdf / cdf[-1]


def histogram_matching(image, params):
    gray = to_gray(image)
    source_counts = np.bincount(gray.ravel(), minlength=256).astype(np.float64)
    source_cdf = np.cumsum(source_counts)
    source_cdf /= source_cdf[-1]
    target_cdf = target_cdf_from_params(params)
    target_levels = np.arange(256)
    mapping = np.interp(source_cdf, target_cdf, target_levels).round().clip(0, 255).astype(np.uint8)
    return gray_to_rgb(mapping[gray])


def local_equalization(image, params):
    size = require_odd(params.get("size", 7), "size")
    gray = to_gray(image)
    radius = size // 2
    padded = np.pad(gray, radius, mode="edge")
    windows = sliding_window_view(padded, (size, size))
    centers = gray[..., None, None]
    ranks = np.count_nonzero(windows <= centers, axis=(-1, -2))
    equalized = np.round(ranks * 255.0 / (size * size)).astype(np.uint8)
    return gray_to_rgb(equalized)


def apply_kernel(image, kernel, normalize_derivative=False):
    source = image.astype(np.float32)
    if image.ndim == 2:
        filtered = cv2.filter2D(source, cv2.CV_32F, kernel, borderType=cv2.BORDER_REFLECT)
    else:
        channels = [cv2.filter2D(source[:, :, idx], cv2.CV_32F, kernel, borderType=cv2.BORDER_REFLECT) for idx in range(source.shape[2])]
        filtered = np.stack(channels, axis=2)
    if normalize_derivative:
        return normalize_to_uint8(filtered)
    return ensure_uint8(np.round(filtered))


def filter_float(image, kernel):
    source = image.astype(np.float32)
    if image.ndim == 2:
        return cv2.filter2D(source, cv2.CV_32F, kernel, borderType=cv2.BORDER_REFLECT)
    channels = [cv2.filter2D(source[:, :, idx], cv2.CV_32F, kernel, borderType=cv2.BORDER_REFLECT) for idx in range(source.shape[2])]
    return np.stack(channels, axis=2)


def box_filter(image, params):
    size = require_odd(params.get("size", 3), "size")
    return cv2.blur(image, (size, size), borderType=cv2.BORDER_REFLECT)


def weighted_average(image, params):
    size = require_odd(params.get("size", 3), "size")
    if "kernel" in params:
        kernel = np.array(params["kernel"], dtype=np.float32)
        if kernel.shape != (size, size):
            raise ProcessingError("kernel dimensions must match size.")
    elif size == 3:
        kernel = np.array([[1, 2, 1], [2, 4, 2], [1, 2, 1]], dtype=np.float32)
    else:
        sigma = max(size / 6.0, 0.1)
        ax = np.arange(-(size // 2), size // 2 + 1, dtype=np.float32)
        xx, yy = np.meshgrid(ax, ax)
        kernel = np.exp(-(xx**2 + yy**2) / (2.0 * sigma**2))
    total = float(np.sum(kernel))
    if math.isclose(total, 0.0):
        raise ProcessingError("weighted average kernel sum must not be zero.")
    return apply_kernel(image, kernel / total)


def median_filter(image, params):
    size = require_odd(params.get("size", 3), "size")
    return cv2.medianBlur(image, size)


def laplacian(image, params):
    mode = params.get("mode", "sharpen")
    lap = filter_float(image, LAPLACIAN_MASK)
    if mode == "edge":
        return normalize_to_uint8(lap)
    sign = params.get("sign", "add")
    source = image.astype(np.float32)
    sharpened = source + lap if sign == "add" else source - lap
    return ensure_uint8(np.round(sharpened))


def high_boost(image, params):
    amplification = float(params.get("amplification", 1.5))
    if not np.isfinite(amplification) or amplification < 1.0:
        raise ProcessingError("amplification must be >= 1.")
    size = require_odd(params.get("size", 3), "size")
    blurred = cv2.blur(image, (size, size), borderType=cv2.BORDER_REFLECT).astype(np.float32)
    boosted = amplification * image.astype(np.float32) - blurred
    return ensure_uint8(np.round(boosted))


def gradient_magnitude(image, gx_kernel, gy_kernel):
    gray = to_gray(image).astype(np.float32)
    gx = cv2.filter2D(gray, cv2.CV_32F, gx_kernel, borderType=cv2.BORDER_REFLECT)
    gy = cv2.filter2D(gray, cv2.CV_32F, gy_kernel, borderType=cv2.BORDER_REFLECT)
    magnitude = np.sqrt(gx**2 + gy**2)
    return gray_to_rgb(normalize_to_uint8(magnitude))


def sobel(image, params):
    return gradient_magnitude(image, SOBEL_GX, SOBEL_GY)


def roberts(image, params):
    return gradient_magnitude(image, ROBERTS_GX, ROBERTS_GY)


def rgb_to_hsi(image):
    rgb = image.astype(np.float32) / 255.0
    r, g, b = rgb[..., 0], rgb[..., 1], rgb[..., 2]
    numerator = 0.5 * ((r - g) + (r - b))
    denominator = np.sqrt((r - g) ** 2 + (r - b) * (g - b)) + 1e-8
    theta = np.arccos(np.clip(numerator / denominator, -1.0, 1.0))
    h = np.where(b <= g, theta, 2.0 * np.pi - theta) / (2.0 * np.pi)
    total = r + g + b
    s = np.where(total <= 1e-8, 0.0, 1.0 - 3.0 * np.minimum(np.minimum(r, g), b) / total)
    i = total / 3.0
    return np.stack([h, s, i], axis=-1)


def hsi_to_rgb(hsi):
    h = (hsi[..., 0] % 1.0) * 2.0 * np.pi
    s = np.clip(hsi[..., 1], 0.0, 1.0)
    i = np.clip(hsi[..., 2], 0.0, 1.0)
    r = np.zeros_like(h)
    g = np.zeros_like(h)
    b = np.zeros_like(h)

    sector0 = h < 2.0 * np.pi / 3.0
    sector1 = (h >= 2.0 * np.pi / 3.0) & (h < 4.0 * np.pi / 3.0)
    sector2 = ~sector0 & ~sector1

    h0 = h[sector0]
    b[sector0] = i[sector0] * (1.0 - s[sector0])
    r[sector0] = i[sector0] * (1.0 + s[sector0] * np.cos(h0) / (np.cos(np.pi / 3.0 - h0) + 1e-8))
    g[sector0] = 3.0 * i[sector0] - (r[sector0] + b[sector0])

    h1 = h[sector1] - 2.0 * np.pi / 3.0
    r[sector1] = i[sector1] * (1.0 - s[sector1])
    g[sector1] = i[sector1] * (1.0 + s[sector1] * np.cos(h1) / (np.cos(np.pi / 3.0 - h1) + 1e-8))
    b[sector1] = 3.0 * i[sector1] - (r[sector1] + g[sector1])

    h2 = h[sector2] - 4.0 * np.pi / 3.0
    g[sector2] = i[sector2] * (1.0 - s[sector2])
    b[sector2] = i[sector2] * (1.0 + s[sector2] * np.cos(h2) / (np.cos(np.pi / 3.0 - h2) + 1e-8))
    r[sector2] = 3.0 * i[sector2] - (g[sector2] + b[sector2])

    return ensure_uint8(np.round(np.clip(np.stack([r, g, b], axis=-1), 0.0, 1.0) * 255.0))


def hsi_intensity_filter(image, params):
    method = params.get("method", "smooth")
    hsi = rgb_to_hsi(image)
    intensity = np.round(hsi[..., 2] * 255.0).astype(np.uint8)
    if method == "sharpen":
        filtered = laplacian(intensity, {"mode": "sharpen", "sign": params.get("sign", "add")})
    else:
        filtered = box_filter(intensity, {"size": params.get("size", 3)})
    hsi[..., 2] = filtered.astype(np.float32) / 255.0
    return hsi_to_rgb(hsi)


def pseudo_color_slices(image, params):
    gray = to_gray(image)
    slices = params.get(
        "slices",
        [
            {"start": 0, "end": 85, "color": [59, 130, 246]},
            {"start": 86, "end": 170, "color": [34, 197, 94]},
            {"start": 171, "end": 255, "color": [239, 68, 68]},
        ],
    )
    output = np.zeros((*gray.shape, 3), dtype=np.uint8)
    for item in slices:
        start = int(item["start"])
        end = int(item["end"])
        color = np.array(item["color"], dtype=np.uint8)
        if start < 0 or end > 255 or start > end or color.shape != (3,):
            raise ProcessingError("Each pseudo-color slice requires start/end in 0..255 and an RGB color.")
        output[(gray >= start) & (gray <= end)] = color
    return output


def gray_to_color_sinusoidal(image, params):
    gray = to_gray(image).astype(np.float32) / 255.0
    hue_frequency = float(params.get("hue_frequency", 1.0))
    saturation_frequency = float(params.get("saturation_frequency", 0.5))
    intensity_frequency = float(params.get("intensity_frequency", 0.25))
    h = (0.5 + 0.5 * np.sin(2.0 * np.pi * hue_frequency * gray)) % 1.0
    s = 0.55 + 0.4 * np.sin(2.0 * np.pi * saturation_frequency * gray + np.pi / 3.0)
    i = 0.5 + 0.45 * np.sin(2.0 * np.pi * intensity_frequency * gray - np.pi / 2.0)
    return hsi_to_rgb(np.stack([h, np.clip(s, 0, 1), np.clip(i, 0, 1)], axis=-1))


OPERATIONS = {
    "negative": negative,
    "log": logarithmic,
    "gamma": gamma,
    "contrast_stretch": contrast_stretch,
    "gray_slice": gray_slice,
    "bit_plane": bit_plane,
    "hist_equalization": histogram_equalization,
    "hist_match": histogram_matching,
    "local_equalization": local_equalization,
    "box_filter": box_filter,
    "weighted_average": weighted_average,
    "median_filter": median_filter,
    "laplacian": laplacian,
    "high_boost": high_boost,
    "sobel": sobel,
    "roberts": roberts,
    "hsi_intensity_filter": hsi_intensity_filter,
    "pseudo_color_slices": pseudo_color_slices,
    "gray_to_color_sinusoidal": gray_to_color_sinusoidal,
}


def process_image(image, operation, params=None):
    params = params or {}
    if operation not in OPERATIONS:
        raise ProcessingError(f"Unsupported operation '{operation}'.")
    return ensure_uint8(OPERATIONS[operation](ensure_uint8(image), params))


def subtract_images(left, right):
    verify_registration([left, right])
    diff = left.astype(np.float32) - right.astype(np.float32)
    return normalize_to_uint8(np.abs(diff))


def average_images(images):
    verify_registration(images)
    stack = np.stack([image.astype(np.float32) for image in images], axis=0)
    return ensure_uint8(np.round(np.mean(stack, axis=0)))


def verify_registration(images):
    if len(images) < 2:
        raise ProcessingError("At least two registered images are required.")
    shape = images[0].shape
    if any(image.shape != shape for image in images[1:]):
        raise ProcessingError("Images must have identical width, height, and channel count.")