import sys
from pathlib import Path

import numpy as np
import pytest
from scipy.sparse import coo_matrix
from scipy.sparse.csgraph import connected_components
from scipy.spatial import cKDTree

sys.path.insert(0, str(Path(__file__).parents[1] / "src"))
from anyimage.operators.cutout_tool import geometry as g
from anyimage.operators.cutout_tool import fine_mesh as fine


def branch_alpha(offset=0):
    yy, xx = np.mgrid[:192, :192]
    xx = xx - offset
    yy = yy - offset
    alpha = (xx - 46) ** 2 + (yy - 48) ** 2 < 30**2
    path = np.asarray(((48, 48), (91, 70), (130, 110), (149, 167)))
    for a, b in zip(path[:-1], path[1:]):
        d = np.stack((xx - a[0], yy - a[1]), axis=-1)
        v = b - a
        t = np.clip(d @ v / (v @ v), 0, 1)
        alpha |= np.sum((d - t[..., None] * v) ** 2, axis=-1) < 2.5**2
    alpha |= (xx - 130) ** 2 + (yy - 110) ** 2 < 7**2
    return alpha.astype(float)


def assert_supported(points, faces, boundary, components):
    triangles = points[faces]
    area = np.abs(
        np.cross(triangles[:, 1] - triangles[:, 0], triangles[:, 2] - triangles[:, 0])
    )
    assert np.isfinite(points).all() and np.all(area > 1e-8)
    assert not boundary[faces].all(axis=1).any()
    edges = np.sort(np.concatenate((faces[:, :2], faces[:, 1:], faces[:, ::2])), axis=1)
    edges, counts = np.unique(edges, axis=0, return_counts=True)
    assert counts.max() <= 2
    assert not boundary[edges[counts == 2]].all(axis=1).any()
    interior = np.flatnonzero(~boundary)
    links = edges[~boundary[edges].any(axis=1)]
    graph = coo_matrix(
        (np.ones(len(links)), (links[:, 0], links[:, 1])),
        shape=(len(points), len(points)),
    ).tocsr()
    count, _labels = connected_components(graph[interior][:, interior], directed=False)
    assert count == len(components)
    profile = g.poisson_balloon_profile(points, faces, boundary)
    assert np.isfinite(profile).all()
    assert np.all(profile[interior] > 0)


@pytest.mark.parametrize("spacing", (4, 8, 16, 28, 32, 36))
@pytest.mark.parametrize("offset", (0, 1))
def test_fine_bent_leg_and_joint_have_connected_positive_support(spacing, offset):
    points, faces, boundary, components = fine.build_fine_mesh(
        branch_alpha(offset), spacing, 0.9
    )
    assert_supported(points, faces, boundary, components)
    interior = points[~boundary]
    for checkpoint in ((90, 70), (130, 110), (146, 159)):
        assert cKDTree(interior).query(np.asarray(checkpoint) + offset)[0] < 12


@pytest.mark.parametrize("spacing", (4, 8, 16, 32))
def test_fine_comb_hole_and_one_pixel_bridge(spacing):
    alpha = np.zeros((160, 180))
    alpha[15:145, 15:42] = 1
    alpha[45:115, 125:165] = 1
    alpha[75, 40:126] = 1
    alpha[65:90, 140:150] = 0
    for row in range(20, 140, 12):
        alpha[row : row + 2, 35:110] = 1
    result = fine.build_fine_mesh(alpha, spacing, 0.9)
    assert_supported(*result)
    points, faces, _boundary, components = result
    assert len(components) == 1
    assert not np.any(
        g.points_inside_contours(
            points[faces].mean(axis=1),
            [np.asarray(((140, 65), (150, 65), (150, 90), (140, 90)))],
        )
    )


def test_fine_cleanup_uses_threshold_and_preserves_input_and_largest_island():
    alpha = np.zeros((40, 60))
    alpha[5:25, 5:25] = 0.5
    alpha[12, 20:45] = 0.5
    alpha[32:34, 50:52] = 1
    original = alpha.copy()
    filtered, stats = fine.filter_alpha(alpha, 0.3)
    assert stats["components_after"] == 1
    assert filtered[12, 44] == 0.5 and filtered[32, 50] == 0
    np.testing.assert_array_equal(alpha, original)
    assert_supported(*fine.build_fine_mesh(alpha, 32, 0.3))
    assert_supported(*fine.build_fine_mesh(alpha, 32, 0.9))
    with pytest.raises(ValueError, match="no visible"):
        fine.build_fine_mesh(np.zeros((10, 10)), 32, 0.9)


def test_higher_detail_recovers_curve_instead_of_splitting_low_edges():
    yy, xx = np.mgrid[:240, :240]
    alpha = ((xx - 120) ** 2 + (yy - 120) ** 2 < 92**2).astype(float)
    reference = g.extract_selection_contours(alpha, max_resolution=240, threshold=0.9)
    source = g.resample_closed_contour(np.asarray(reference.components[0][0]), 0.5)
    errors = []
    for spacing in (32, 8, 4):
        contours = fine.contours_for_component(reference.components[0], spacing)
        sampled = g.resample_closed_contour(contours[0], 0.5)
        forward = fine.ContourDomain(contours).clearance(source)
        backward = fine.ContourDomain(
            [np.asarray(reference.components[0][0])]
        ).clearance(sampled)
        errors.append(max(np.percentile(forward, 95), np.percentile(backward, 95)))
    assert errors[1] < errors[0]
    assert errors[2] <= errors[1] + 0.1


def test_segment_validation_rejects_holes_tangency_and_outside():
    domain = fine.ContourDomain(
        [
            np.asarray(((0, 0), (20, 0), (20, 20), (0, 20))),
            np.asarray(((8, 8), (12, 8), (12, 12), (8, 12))),
        ]
    )
    starts = np.asarray(((2, 2), (2, 10), (2, 8), (-1, 2), (2, 2)))
    ends = np.asarray(((18, 2), (18, 10), (18, 8), (2, 2), (2, 2)))
    assert domain.valid_segments(starts, ends).tolist() == [
        True,
        False,
        False,
        False,
        True,
    ]


def test_closed_mode_uses_original_geometry_without_fine_import(monkeypatch):
    from types import SimpleNamespace
    from mathutils import Matrix
    from anyimage.common.image import image_empty_bounds

    source = SimpleNamespace(
        data=SimpleNamespace(size=(80, 60)),
        empty_display_size=2.0,
        empty_image_offset=(-0.5, -0.5),
        matrix_world=Matrix.Identity(4),
    )
    monkeypatch.setattr(
        g,
        "addon_preferences",
        lambda _context: SimpleNamespace(balloon_profile_algorithm="POISSON"),
    )
    monkeypatch.setattr(
        fine,
        "build_fine_mesh",
        lambda *_args: pytest.fail("Fine algorithm called while disabled"),
    )
    alpha = np.ones((60, 80))
    actual = g.build_base_shape(
        None, source, 16, alpha, (0, 0, 80, 60), fine_outline=False
    )
    contours = g.extract_cutout_contours(alpha, 16, threshold=0.9)
    expected = g.cutout_domain(
        contours,
        image_empty_bounds(source),
        g.cutout_sampling_spacing(source, 16),
        ((1.0, 0.0), (0.0, 1.0)),
    )
    np.testing.assert_array_equal(actual.vertices, expected[0])
    np.testing.assert_array_equal(actual.uv, expected[2])
    np.testing.assert_array_equal(actual.balloon, expected[-1])


def test_width_index_matches_exhaustive_covering_disks():
    from scipy import ndimage

    width = fine.WidthSamples(branch_alpha().astype(bool))
    points = np.asarray(((48, 48), (85, 70), (130, 110), (147, 160), (2, 2)))
    radii = []
    for point in points:
        covered = np.sum((width.centers - point) ** 2, axis=1) <= width.reach_squared
        radii.append(max(0.5, width.radius[covered].max() if covered.any() else 0.5))
    distances = ndimage.map_coordinates(
        width.distance,
        [points[:, 1] - 0.5, points[:, 0] - 0.5],
        order=1,
        mode="nearest",
    )
    expected = np.maximum(
        0.7, np.minimum(32, np.minimum(16, np.asarray(radii) * 0.9) + 0.4 * distances)
    )
    np.testing.assert_allclose(width(points), expected)


def test_relaxation_reduces_sampling_error_without_adding_points():
    contours = [np.asarray(((4.0, 4.0), (194.0, 4.0), (194.0, 104.0), (4.0, 104.0)))]
    spine = np.empty((0, 2))
    original = g.uniform_interior_points_for_contours(contours, 16)
    relaxed = fine.relax_points(contours, spine, (110, 200), 16)
    yy, xx = np.mgrid[6:104:4, 6:194:4]
    samples = np.column_stack((xx.ravel() + 0.5, yy.ravel() + 0.5))

    def energy(free):
        return np.mean(
            cKDTree(np.concatenate((*contours, free))).query(samples)[0] ** 2
        )

    assert len(relaxed) == len(original)
    assert energy(relaxed) < energy(original)
    assert fine.ContourDomain(contours).contains(relaxed).all()


def test_cdt_preserves_constraint_split_by_an_interior_vertex():
    contours = [np.asarray(((0.0, 0.0), (20.0, 0.0), (20.0, 10.0), (0.0, 10.0)))]
    spine = np.asarray(((2.0, 5.0), (18.0, 5.0)))
    points, faces, constraints = fine.triangulate(
        contours, spine, [(0, 1)], np.asarray(((10.0, 5.0),)), 8
    )
    assert len(constraints) == 2
    assert np.isclose(
        np.linalg.norm(
            points[constraints[:, 1]] - points[constraints[:, 0]], axis=1
        ).sum(),
        16,
    )
    assert not g.boundary_vertices_from_faces(len(points), faces)[constraints].any()


def test_contour_validation_accepts_pixel_contacts_and_rejects_crossings():
    crossing = np.asarray(((0, 0), (4, 4), (0, 4), (4, 0)))
    assert not fine.valid_contours([crossing])
    touching = np.asarray(
        ((0, 0), (2, 0), (2, 2), (4, 2), (4, 4), (2, 4), (2, 2), (0, 2))
    )
    assert fine.valid_contours([touching])


@pytest.mark.parametrize("offset", [0, 2048, 4096, 8192, 16384])
@pytest.mark.parametrize("spacing", [4, 0.5])
def test_cdt_preserves_short_constraints_at_large_pixel_coordinates(offset, spacing):
    contour = np.asarray(((0.0, 0.0), (20.0, 0.0), (20.0, 20.0), (0.0, 20.0))) + offset
    spine = np.asarray(((5.1234567, 5.654321), (6.456789, 6.234567))) + offset
    free = g.uniform_interior_points_for_contours([contour], spacing)
    points, faces, constraints = fine.triangulate(
        [contour], spine, [(0, 1)], free, spacing
    )
    assert len(constraints) >= 1
    assert not g.boundary_vertices_from_faces(len(points), faces)[constraints].any()


@pytest.mark.parametrize(
    "damage", ["missing_edge", "missing_faces", "missing_endpoint"]
)
def test_cdt_rejects_broken_constraint_coverage(monkeypatch, damage):
    real_cdt = fine.delaunay_2d_cdt

    def damaged_cdt(*args):
        result = list(real_cdt(*args))
        result[2] = list(result[2])
        result[3] = [list(ids) for ids in result[3]]
        result[4] = [list(ids) for ids in result[4]]
        constraint = len(args[1]) - 1
        edge_id = next(i for i, ids in enumerate(result[4]) if constraint in ids)
        if damage == "missing_edge":
            result[4][edge_id].remove(constraint)
        elif damage == "missing_faces":
            a, b = result[1][edge_id]
            result[2] = [face for face in result[2] if not (a in face and b in face)]
        else:
            endpoint = args[1][-1][1]
            for originals in result[3]:
                if endpoint in originals:
                    originals.remove(endpoint)
        return result

    monkeypatch.setattr(fine, "delaunay_2d_cdt", damaged_cdt)
    contour = np.asarray(((0.0, 0.0), (20.0, 0.0), (20.0, 10.0), (0.0, 10.0)))
    spine = np.asarray(((2.0, 5.0), (18.0, 5.0)))
    with pytest.raises(ValueError, match="lost an interior constraint"):
        fine.triangulate([contour], spine, [(0, 1)], np.asarray(((10.0, 5.0),)), 8)