remesh 0.0.5

// SPDX-License-Identifier: MIT OR Apache-2.0
// Copyright (c) 2025 lacklustr@protonmail.com https://github.com/eadf

use crate::common::VertexIndex;
use crate::corner_table::{CornerIndex, TriangleIndex};
use crate::isotropic_remesh::IsotropicRemeshAlgo;
use crate::prelude::*;
use vector_traits::glam::{Mat3, Vec3, Vec3A};

#[test]
fn test_coplanar_triangles_same_orientation() {
    // Two triangles in the XY plane, both CCW when viewed from +Z
    // Triangle 1 (abc): (0,0,0) -> (1,0,0) -> (0,1,0)
    // Triangle 2 (dba): (0,-1,0) -> (1,0,0) -> (0,0,0)  [shared edge (1,0,0)->(0,0,0)]
    let a = Vec3A::new(0.0, 0.0, 0.0);
    let b = Vec3A::new(1.0, 0.0, 0.0);
    let c = Vec3A::new(0.0, 1.0, 0.0);
    let d = Vec3A::new(0.0, -1.0, 0.0);

    let coplanarity =
        IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_opposite_winding_slow_wrong_order(
            a, b, c, d,
        );

    println!("Coplanar same orientation: {:.6}", coplanarity);

    // Should be close to 1.0 (parallel normals)
    assert!(
        (coplanarity - 1.0).abs() < 0.001,
        "Coplanar triangles should have coplanarity ≈ 1.0 (coplanarity={coplanarity})"
    );
    let coplanarity =
        IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_same_winding(a, b, c, d, 0.95 * 0.95);
    assert_eq!(coplanarity, Some(true));
}

#[test]
fn test_45_degree_dihedral_angle() {
    // Two triangles forming a 45-degree dihedral angle
    // Triangle 1 (ABC): in XY plane, normal pointing in +Z
    //   A(0,0,0) -> B(1,0,0) -> C(0,1,0)
    //   Normal: (0, 0, 1)
    // Triangle 2 (DBA): tilted 45° from XY plane
    //   D rotated 45° around the shared edge (X-axis)
    //   D(0, -1/√2, 1/√2) creates a 45° angle
    let a = Vec3A::new(0.0, 1.0, 0.0);
    let b = Vec3A::new(0.0, 0.0, 0.0);
    let c = Vec3A::new(1.0, 0.0, 0.0);
    let sqrt2_inv = std::f32::consts::FRAC_1_SQRT_2;
    let d = Vec3A::new(0.0, -sqrt2_inv, sqrt2_inv);

    let coplanarity =
        IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_opposite_winding_slow(a, b, c, d);

    println!("45-degree dihedral angle: {:.6}", coplanarity);

    // Should be close to cos(45°) ≈ 0.707
    assert!(
        (coplanarity - sqrt2_inv).abs() < 0.01,
        "45° dihedral should have coplanarity ≈ 0.707 (coplanarity={coplanarity} ({}°))",
        coplanarity.acos().to_degrees()
    );

    // With 45-degree threshold - should be right at the boundary
    let threshold_45_deg_sq = sqrt2_inv * sqrt2_inv;
    let result = IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_opposite_winding(
        a,
        b,
        c,
        d,
        threshold_45_deg_sq,
    );

    assert_eq!(
        result,
        Some(true),
        "45° dihedral with 45° threshold should pass (at boundary)"
    );
}

#[test]
fn test_90_degree_triangle_angle() {
    // Two triangles forming a 90-degree dihedral angle
    let a = Vec3A::new(0.0, 0.0, 1.0);
    let b = Vec3A::new(0.0, 0.0, 0.0);
    let c = Vec3A::new(0.0, 1.0, 0.0);
    let d = Vec3A::new(1.0, 0.0, 0.0);
    let sqrt2_inv = std::f32::consts::FRAC_1_SQRT_2;

    let coplanarity =
        IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_opposite_winding_slow_wrong_order(
            a, b, c, d,
        );

    println!("45-degree dihedral angle: {:.6}", coplanarity);

    // Should be close to cos(45°) ≈ 0.0
    assert!(
        (coplanarity).abs() < 0.01,
        "90° dihedral should have coplanarity ≈ 0.0 (coplanarity={coplanarity})"
    );

    // With 45-degree threshold - 90° triangles should be well outside the boundary
    let threshold_45_deg_sq = sqrt2_inv * sqrt2_inv;
    let result = IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_opposite_winding(
        a,
        b,
        c,
        d,
        threshold_45_deg_sq,
    );

    // this assert fails! result == Some(true)
    assert_eq!(
        result, None,
        "45° dihedral with 45° threshold should *not* pass for 90° triangles"
    );
}

#[test]
fn test_dihedral_angle_rotation() {
    // Base vertices
    let b = Vec3A::new(0.0, 0.0, 0.0);
    let c = Vec3A::new(0.0, 1.0, 0.0);
    let d = Vec3A::new(1.0, 0.0, 0.0);

    // Test angles from 90 degrees down to 2 degrees in 5-degree increments
    for degrees in (2..=90).rev().step_by(5) {
        let radians = (degrees as f32).to_radians();

        // Rotate point 'a' around Y-axis from Z-axis toward X-axis
        // At 0°: a = (0, 0, 1) - perpendicular to triangle bcd
        // At 90°: a = (1, 0, 0) - coplanar with triangle bcd
        let a = Vec3A::new(radians.cos(), 0.0, radians.sin());

        let coplanarity =
            IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_same_winding_slow(a, b, c, d);
        let expected_cos = radians.cos(); // Theoretical cosine of the dihedral angle

        println!(
            "Testing:{}° dihedral angle: cos_actual={:.6}, cos_expected={:.6}",
            degrees, coplanarity, expected_cos
        );
        assert!(
            (coplanarity - expected_cos).abs() < 0.001,
            "expected_cos!=coplanarity {expected_cos}!={coplanarity} coplanarity={}°",
            coplanarity.acos().to_degrees()
        );

        // Test with threshold slightly above and below the actual angle
        let threshold_above = (radians + 1.0f32.to_radians()).cos().powi(2);
        let threshold_below = (radians - 1.0f32.to_radians()).cos().powi(2);
        let result_above = IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_opposite_winding(
            a,
            b,
            c,
            d,
            threshold_above,
        );
        assert_ne!(
            result_above,
            Some(true),
            "dihedral angle threshold {}° should accept triangles set at {degrees}°",
            threshold_above.sqrt().acos().to_degrees()
        );
        let result_below = IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_opposite_winding(
            a,
            b,
            c,
            d,
            threshold_below,
        );
        assert_eq!(
            result_below,
            None,
            "dihedral angle threshold {}° should reject triangles set at {degrees}°",
            threshold_below.sqrt().acos().to_degrees()
        );
        //println!("  Threshold +1°: {:?}, -1°: {:?}", result_above, result_below);
    }
}

#[test]
fn test_dihedral_angle_rotation_2() {
    // Base vertices
    let b = Vec3A::new(0.0, 0.0, 0.0);
    let c = Vec3A::new(0.0, 1.0, 0.0);
    let d = Vec3A::new(1.0, 0.0, 0.0);

    // Test angles from 90° down to 2° in 5-degree increments
    for degrees in (-90_i32..=90).step_by(2) {
        let radians = (degrees as f32).to_radians();

        // Rotate point 'a' around Y-axis from Z-axis toward X-axis
        // At 0°: a = (0, 0, 1) - perpendicular to triangle bcd
        // At 90°: a = (-1, 0, 0) - coplanar with triangle bcd
        let a = Vec3A::new(-radians.cos(), 0.0, radians.sin());

        let coplanarity =
            IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_opposite_winding_slow(a, b, c, d);
        let expected_cos = radians.cos(); // Theoretical cosine of the dihedral angle

        println!(
            "Testing:{}° dihedral angle: cos_actual={:.6}, cos_expected={:.6} a:{:?}",
            degrees,
            coplanarity,
            expected_cos,
            <Vec3A as Into<[f32; 3]>>::into(a)
        );
        assert!(
            (coplanarity - expected_cos).abs() < 0.001,
            "expected_cos!=coplanarity {expected_cos}!={coplanarity} coplanarity={}°",
            coplanarity.acos().to_degrees()
        );

        let threshold_sq_above = (radians.abs() - 1.0f32.to_radians()).cos().powi(2);
        let threshold_sq_below = (radians.abs() + 1.0f32.to_radians()).cos().powi(2);
        println!(
            "{degrees}°.abs()={}° above:{}° below:{}° coplanarity^2:{} threshold_sq_above:{threshold_sq_above} threshold_sq_below:{threshold_sq_below}",
            degrees.abs(),
            threshold_sq_above.sqrt().acos().to_degrees(),
            threshold_sq_below.sqrt().acos().to_degrees(),
            coplanarity.powi(2)
        );
    }
}

#[test]
fn test_dihedral_angle_thresholds() {
    let b = Vec3A::new(0.0, 0.0, 0.0);
    let c = Vec3A::new(0.0, 1.0, 0.0);
    let d = Vec3A::new(1.0, 0.0, 0.0);

    // Test various threshold angles
    for threshold_degrees in (3..=85).step_by(2) {
        let threshold_radians = (threshold_degrees as f32).to_radians();
        let threshold_cos_sq = threshold_radians.cos().powi(2);

        // Create test cases that should barely PASS (angle slightly smaller than threshold)
        let pass_angle = threshold_degrees - 1;
        let pass_radians = (pass_angle as f32).to_radians();
        let a_pass = Vec3A::new(-pass_radians.cos(), 0.0, pass_radians.sin());
        let cos_sq_pass =
            IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_opposite_winding_slow(
                a_pass, b, c, d,
            )
            .powi(2);

        // Should PASS: cos_sq > threshold_cos_sq (angle < threshold)
        assert!(
            cos_sq_pass > threshold_cos_sq,
            "Pass case failed: {}° should pass {}° threshold",
            pass_angle,
            threshold_degrees
        );

        // Create test cases that should barely FAIL (angle slightly larger than threshold)
        let fail_angle = threshold_degrees + 1;
        let fail_radians = (fail_angle as f32).to_radians();
        let a_fail = Vec3A::new(-fail_radians.cos(), 0.0, fail_radians.sin());
        let cos_sq_fail =
            IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_opposite_winding_slow(
                a_fail, b, c, d,
            )
            .powi(2);

        // Should FAIL: cos_sq < threshold_cos_sq (angle > threshold)
        assert!(
            cos_sq_fail < threshold_cos_sq,
            "Fail case failed: {}° should fail {}° threshold",
            fail_angle,
            threshold_degrees
        );

        println!(
            "Threshold {}°: PASS at {}° (cos_sq={:.6}), FAIL at {}° (cos_sq={:.6})",
            threshold_degrees, pass_angle, cos_sq_pass, fail_angle, cos_sq_fail
        );
    }
}

#[test]
fn test_dihedral_angle_thresholds_works() {
    let b = Vec3A::new(0.0, 0.0, 0.0);
    let c = Vec3A::new(0.0, 1.0, 0.0);
    let d = Vec3A::new(1.0, 0.0, 0.0);

    for threshold_degrees in (2..=88).step_by(1) {
        let threshold_radians = (threshold_degrees as f32).to_radians();
        let threshold_cos_sq = threshold_radians.cos().powi(2);

        // Test positive dihedral angles
        let pass_angle = threshold_degrees - 1;
        let fail_angle = threshold_degrees + 1;

        let pass_radians = (pass_angle as f32).to_radians();
        let fail_radians = (fail_angle as f32).to_radians();

        // Positive dihedral (a on positive Z side)
        let a_pass_pos = Vec3A::new(-pass_radians.cos(), 0.0, pass_radians.sin());
        let a_fail_pos = Vec3A::new(-fail_radians.cos(), 0.0, fail_radians.sin());

        let cos_sq_pass_pos =
            IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_opposite_winding_slow(
                a_pass_pos, b, c, d,
            )
            .powi(2);
        let cos_sq_fail_pos =
            IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_opposite_winding_slow(
                a_fail_pos, b, c, d,
            )
            .powi(2);

        assert!(
            cos_sq_pass_pos > threshold_cos_sq,
            "Positive pass case failed"
        );
        assert!(
            cos_sq_fail_pos < threshold_cos_sq,
            "Positive fail case failed"
        );

        // Negative dihedral (a on negative Z side) - should give same results
        let a_pass_neg = Vec3A::new(-pass_radians.cos(), 0.0, -pass_radians.sin());
        let a_fail_neg = Vec3A::new(-fail_radians.cos(), 0.0, -fail_radians.sin());

        assert_eq!(
            Some(true),
            IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_opposite_winding(
                a_pass_neg,
                b,
                c,
                d,
                threshold_cos_sq
            )
        );
        assert_eq!(
            Some(true),
            IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_opposite_winding(
                a_pass_pos,
                b,
                c,
                d,
                threshold_cos_sq
            )
        );
        assert_eq!(
            None,
            IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_opposite_winding(
                a_fail_neg,
                b,
                c,
                d,
                threshold_cos_sq
            )
        );
        assert_eq!(
            None,
            IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_opposite_winding(
                a_fail_pos,
                b,
                c,
                d,
                threshold_cos_sq
            )
        );

        println!(
            "Threshold {}°: PASS both signs at {}° (cos_sq={:.6}), FAIL both signs at {}° (cos_sq={:.6})",
            threshold_degrees, pass_angle, cos_sq_pass_pos, fail_angle, cos_sq_fail_pos
        );
    }
}

#[test]
fn test_opposite_normals() {
    // Two triangles where one is flipped relative to the other
    // Triangle 1 (abc): (0,0,0) -> (1,0,0) -> (0,1,0)  [normal points +Z]
    // Triangle 2 (dba): (0,2,0) -> (1,0,0) -> (0,0,0)  [normal points -Z]
    let a = Vec3A::new(0.0, 0.0, 0.0);
    let b = Vec3A::new(1.0, 0.0, 0.0);
    let c = Vec3A::new(0.0, 1.0, 0.0);
    let d = Vec3A::new(0.0, 2.0, 0.0);

    let coplanarity =
        IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_opposite_winding_slow_wrong_order(
            a, b, c, d,
        );

    println!("Opposite normals: {:.6}", coplanarity);

    let abc = (b - a).cross(c - a);
    let dba = (b - d).cross(a - d);

    println!("  Triangle abc normal: {:?}", abc.normalize());
    println!("  Triangle dba normal: {:?}", dba.normalize());

    // Should be close to -1.0 (opposite normals)
    assert!(
        (coplanarity + 1.0).abs() < 0.001,
        "Opposite triangles should have coplanarity ≈ -1.0"
    );
}

#[test]
fn test_perpendicular_triangles() {
    // Two triangles at 90 degrees to each other
    // Triangle 1 (abc): in XY plane -> normal points +Z
    // Triangle 2 (dba): in XZ plane -> normal points +Y
    let a = Vec3A::new(0.0, 0.0, 0.0);
    let b = Vec3A::new(1.0, 0.0, 0.0);
    let c = Vec3A::new(0.0, 1.0, 0.0); // XY plane triangle
    let d = Vec3A::new(0.0, 0.0, 1.0); // Makes XZ plane triangle

    let coplanarity =
        IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_opposite_winding_slow_wrong_order(
            a, b, c, d,
        );

    println!("Perpendicular triangles: {:.6}", coplanarity);

    let normal1 = (b - a).cross(c - a); // (1,0,0) x (0,1,0) = (0,0,1)
    let normal2 = (d - b).cross(a - b); // (-1,0,1) x (-1,0,0) = ?

    println!("  Triangle abc normal: {:?}", normal1.normalize());
    println!("  Triangle dba normal: {:?}", normal2.normalize());

    // Should be close to 0.0 (perpendicular normals)
    assert!(
        coplanarity.abs() < 0.001,
        "Perpendicular triangles should have coplanarity ≈ 0.0"
    );
}

#[test]
fn test_shared_edge_ccw_triangles() {
    // Create two adjacent triangles that share an edge, both CCW

    //     c
    //    /|\
    //   / | \
    //  /  |  \
    // a---+---b
    //  \  |  /
    //   \ | /
    //    \|/
    //     d

    let a = Vec3A::new(0.0, 0.0, 0.0); // shared vertex 1
    let b = Vec3A::new(2.0, 0.0, 0.0); // shared vertex 2
    let c = Vec3A::new(1.0, 1.0, 0.0); // third vertex of triangle 1
    let d = Vec3A::new(1.0, -1.0, 0.0); // third vertex of triangle 2

    // Triangle 1 (abc): a->b->c should be CCW when viewed from +Z
    // Triangle 2 (dba): d->b->a should be CCW when viewed from +Z
    let coplanarity =
        IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_opposite_winding_slow_wrong_order(
            a, b, c, d,
        );

    println!("Shared edge CCW triangles: {:.6}", coplanarity);

    // Calculate normals manually to verify
    let normal1 = (b - a).cross(c - a); // (2,0,0) x (1,1,0) = (0,0,2)
    let normal2 = (b - d).cross(a - d).normalize(); // (2,0,0) x (1,-1,0) = (0,0,2)

    println!(
        "  Triangle abc normal: {:?} -> {:?}",
        normal1,
        normal1.normalize()
    );
    println!(
        "  Triangle dba normal: {:?} -> {:?}",
        normal2,
        normal2.normalize()
    );

    // These should be coplanar but opposite facing (this is the problem case!)
    // If both triangles are truly CCW in a valid mesh, they should have similar normals
    // But the math here gives opposite normals because of how the shared edge is oriented

    println!("  This case shows the issue: both triangles are CCW but normals are opposite");
    assert!(
        (coplanarity - 1.0).abs() < 0.001,
        "This should give 1.0, demonstrating the issue coplanarity:{coplanarity}"
    );
}

#[test]
fn test_proper_adjacent_triangles() {
    // Let's create the correct case for adjacent triangles in a mesh
    // Both triangles share edge a-b, and both have outward-facing normals

    //      c
    //     /|
    //    / |
    //   /  |
    //  a---b
    //   \  |
    //    \ |
    //     \|
    //      d

    let a = Vec3A::new(0.0, 0.0, 0.0);
    let b = Vec3A::new(1.0, 0.0, 0.0);
    let c = Vec3A::new(0.5, 1.0, 0.1); // slightly out of plane
    let d = Vec3A::new(0.5, -1.0, -0.1); // on the opposite side

    // For a proper mesh, if we're looking at the "outside" surface:
    // Triangle 1: a->b->c (CCW from outside)
    // Triangle 2: b->a->d (CCW from outside)
    // This means for the method call triangles_coplanarity(a,b,c,d):
    // - Triangle abc: a->b->c
    // - Triangle dba: d->b->a

    let coplanarity =
        IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_opposite_winding_slow_wrong_order(
            a, b, c, d,
        );

    println!("Proper adjacent triangles: {:.6}", coplanarity);

    let normal1 = (b - a).cross(c - a);
    let normal2 = (d - b).cross(a - b);

    println!(
        "  Triangle abc normal: {:?} -> {:?}",
        normal1,
        normal1.normalize()
    );
    println!(
        "  Triangle dba normal: {:?} -> {:?}",
        normal2,
        normal2.normalize()
    );

    // These should have similar normals (both pointing roughly "outward")
    assert!(
        coplanarity > 0.5,
        "Properly oriented adjacent triangles should have similar normals"
    );
}

#[test]
fn debug_specific_case() {
    let vc = Vec3A::new(-1.3545021, -1.0532312, -0.2397142); // V0
    let vn = Vec3A::new(-0.8545021, -1.0532312, 0.2602858); // V2
    let vp = Vec3A::new(-1.3545021, -1.0532312, 0.2602858); // V5
    let vo = Vec3A::new(-0.8545021, -1.0532312, -0.2397142); // V6

    println!("  vc: {:?}", vc);
    println!("  vn: {:?}", vn);
    println!("  vp: {:?}", vp);
    println!("  vo: {:?}", vo);

    let coplanarity =
        IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_opposite_winding_slow_wrong_order(
            vc, vn, vp, vo,
        );

    println!("  Coplanarity: {:.6}", coplanarity);

    let normal1 = (vn - vc).cross(vp - vc);
    let normal2 = (vo - vn).cross(vc - vn);

    println!("  Normal1: {:?} -> {:?}", normal1, normal1.normalize());
    println!("  Normal2: {:?} -> {:?}", normal2, normal2.normalize());

    // Let's see what happens if we try different vertex orderings
    println!("\n  Trying different orderings:");

    // Option 1: Reverse triangle 2
    let coplanarity2 =
        IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_opposite_winding_slow_wrong_order(
            vn, vc, vp, vo,
        );
    println!(
        "    triangles_coplanarity(vn, vc, vp, vo): {:.6}",
        coplanarity2
    );

    // Option 2: Different arrangement
    let coplanarity3 =
        IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_opposite_winding_slow_wrong_order(
            vc, vn, vp, vo,
        );
    println!(
        "    Original call triangles_coplanarity(vc, vn, vp, vo): {:.6}",
        coplanarity3
    );
}

#[test]
fn test_coplanar_triangles_same_orientatio_2() {
    // Two triangles in the XY plane, both CCW when viewed from +Z
    // Shared edge: (b,c)
    let a = Vec3A::new(0.0, 0.0, 0.0);
    let b = Vec3A::new(0.0, 1.0, 0.0);
    let c = Vec3A::new(1.0, 1.0, 0.0);
    let d = Vec3A::new(1.0, 0.0, 0.0);
    let threshold_sq = 0.9_f32.powi(2);

    let result =
        IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_same_winding(a, b, c, d, threshold_sq);
    assert_eq!(result, Some(true));
}

#[test]
fn test_coplanar_triangles_opposite_orientation() {
    // Same plane, but one flipped (reverse winding)
    let a = Vec3A::new(0.0, 0.0, 0.0);
    let b = Vec3A::new(0.0, 1.0, 0.0);
    let c = Vec3A::new(1.0, 1.0, 0.0);
    let d = Vec3A::new(1.0, 0.0, 0.0);
    let threshold_sq = 0.9_f32.powi(2);

    // Flip winding by swapping order of b/c in one triangle
    let result =
        IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_same_winding(a, c, b, d, threshold_sq);
    assert_eq!(result, Some(true));
}

#[test]
fn test_triangles_tilted_small_angle() {
    // Tilt DBC triangle 10° around shared edge BC
    let angle = 10_f32.to_radians();
    let a = Vec3A::new(0.0, 0.0, 0.0);
    let b = Vec3A::new(0.0, 1.0, 0.0);
    let c = Vec3A::new(1.0, 1.0, 0.0);
    let d = Vec3A::new(1.0, 0.0, angle.sin());
    let threshold_sq = angle.cos().powi(2);

    let result =
        IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_same_winding(a, b, c, d, threshold_sq);
    assert_eq!(result, Some(true));
}

#[test]
fn test_triangles_tilted_large_angle() {
    // Tilt DBC triangle 60° around shared edge BC
    let angle = 60_f32.to_radians();
    let a = Vec3A::new(0.0, 0.0, 0.0);
    let b = Vec3A::new(0.0, 1.0, 0.0);
    let c = Vec3A::new(1.0, 1.0, 0.0);
    let d = Vec3A::new(1.0, 0.0, angle.sin());
    let threshold_sq = (30_f32.to_radians()).cos().powi(2); // accept up to 30°

    let result =
        IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_same_winding(a, b, c, d, threshold_sq);
    assert_eq!(result, None);
}

#[test]
fn test_degenerate_triangle() {
    let a = Vec3A::new(0.0, 1.0, 0.0); // degenerate ABC
    let b = Vec3A::new(0.0, 1.0, 0.0);
    let c = Vec3A::new(1.0, 1.0, 0.0);
    let d = Vec3A::new(0.0, 0.0, 0.0);

    let threshold_sq = 0.9_f32.powi(2);
    let result =
        IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_same_winding(a, b, c, d, threshold_sq);
    assert_eq!(result, Some(false));
}

#[test]
fn test_ok_triangle_2() {
    let threshold = 0.9;

    let a: Vec3A = [0.0, 0.0, 0.0].into();
    let b: Vec3A = [1.0, 0.0, 0.0].into();
    let c: Vec3A = [1.0, 1.0, 0.0].into();
    let d: Vec3A = [0.0, 1.0, 0.0].into();

    assert!(IsotropicRemeshAlgo::<f32, Vec3>::quad_is_convex_coplanar(
        a, b, c, d, threshold
    ));
}

#[test]
fn test_arrowhead_quad_1() {
    let threshold = 0.9;

    let a: Vec3A = [0.0, 0.0, 0.0].into(); // all planar and convex
    let b: Vec3A = [1.0, 0.0, 0.0].into();
    let c: Vec3A = [1.0, 1.0, 0.0].into();
    let d: Vec3A = [0.0, 1.0, 0.0].into();
    assert!(IsotropicRemeshAlgo::<f32, Vec3>::quad_is_convex_coplanar(
        a, b, c, d, threshold
    ));

    let a: Vec3A = [0.9, 0.9, 0.0].into(); // 'a' almost at c, should fail
    let b: Vec3A = [1.0, 0.0, 0.0].into();
    let c: Vec3A = [1.0, 1.0, 0.0].into();
    let d: Vec3A = [0.0, 1.0, 0.0].into();
    assert!(!IsotropicRemeshAlgo::<f32, Vec3>::quad_is_convex_coplanar(
        a, b, c, d, threshold
    ));

    let a: Vec3A = [0.0, 0.0, 0.0].into();
    let b: Vec3A = [0.1, 0.9, 0.0].into(); // b almost at d, should fail
    let c: Vec3A = [1.0, 1.0, 0.0].into();
    let d: Vec3A = [0.0, 1.0, 0.0].into();
    assert!(!IsotropicRemeshAlgo::<f32, Vec3>::quad_is_convex_coplanar(
        a, b, c, d, threshold
    ));

    let a: Vec3A = [0.0, 0.0, 0.0].into();
    let b: Vec3A = [1.0, 0.0, 0.0].into();
    let c: Vec3A = [0.1, 0.1, 0.0].into(); // c almost at a, should fail
    let d: Vec3A = [0.0, 1.0, 0.0].into();
    assert!(!IsotropicRemeshAlgo::<f32, Vec3>::quad_is_convex_coplanar(
        a, b, c, d, threshold
    ));

    let a: Vec3A = [0.0, 0.0, 0.0].into();
    let b: Vec3A = [1.0, 0.0, 0.0].into();
    let c: Vec3A = [1.0, 1.0, 0.0].into();
    let d: Vec3A = [0.9, 0.1, 0.0].into(); // d almost at b, should fail
    assert!(!IsotropicRemeshAlgo::<f32, Vec3>::quad_is_convex_coplanar(
        a, b, c, d, threshold
    ));

    // Test 1: Perfect planar convex quad (should pass)
    let a: Vec3A = [0.0, 0.0, 0.0].into();
    let b: Vec3A = [1.0, 0.0, 0.0].into();
    let c: Vec3A = [1.0, 1.0, 0.0].into();
    let d: Vec3A = [0.0, 1.0, 0.0].into();
    assert!(IsotropicRemeshAlgo::<f32, Vec3>::quad_is_convex_coplanar(
        a, b, c, d, threshold
    ));

    // Test 2: Non-axis-aligned convex planar quad (should pass)
    let a: Vec3A = [0.0, 0.0, 0.0].into();
    let b: Vec3A = [1.0, 0.2, 0.0].into();
    let c: Vec3A = [0.8, 1.0, 0.0].into();
    let d: Vec3A = [-0.2, 0.8, 0.0].into();
    assert!(IsotropicRemeshAlgo::<f32, Vec3>::quad_is_convex_coplanar(
        a, b, c, d, threshold
    ));

    // Test 3: Concave quad (should fail)
    let a: Vec3A = [0.0, 0.0, 0.0].into();
    let b: Vec3A = [1.0, 0.0, 0.0].into();
    let c: Vec3A = [0.2, 0.5, 0.0].into(); // Creates concavity
    let d: Vec3A = [0.0, 1.0, 0.0].into();
    assert!(!IsotropicRemeshAlgo::<f32, Vec3>::quad_is_convex_coplanar(
        a, b, c, d, threshold
    ));

    // Test 5: Highly non-planar quad (should fail)
    let a: Vec3A = [0.0, 0.0, 0.0].into();
    let b: Vec3A = [1.0, 0.0, 0.0].into();
    let c: Vec3A = [1.0, 1.0, 1.0].into(); // Far out of plane
    let d: Vec3A = [0.0, 1.0, 0.0].into();
    assert!(!IsotropicRemeshAlgo::<f32, Vec3>::quad_is_convex_coplanar(
        a, b, c, d, threshold
    ));
}

#[test]
fn dihedral_angle_cosine_test_negative_threshold() {
    let a: Vec3A = [0.0, 0.0, 1.0].into();
    let b: Vec3A = [1.0, 0.0, 0.0].into();
    let c: Vec3A = [0.5, 0.0, 1.0].into();
    let d: Vec3A = [0.5, 0.5, 0.5].into();

    let threshold_sq = (-0.9_f32).powi(2);
    let result =
        IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_crease_angle(a, b, c, d, threshold_sq);
    assert_eq!(result, Some(true));
}

#[test]
fn preserves_surface_normals() {
    let a: Vec3A = [0.5, 0.0, 1.0].into();
    let b: Vec3A = [0.0, 0.0, 1.0].into(); // 55°
    let c: Vec3A = [1.0, 0.0, 0.0].into();
    let d: Vec3A = [0.3, 0.3, 1.0].into();

    let slow_dihedral =
        IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_same_winding_slow(a, b, c, d);

    let threshold_sq = 56.0_f32.to_radians().cos().powi(2);
    let result =
        IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_same_winding(a, b, c, d, threshold_sq);
    assert_eq!(
        result,
        Some(true),
        "failed cos at {}° slow_dihedral:{}°",
        threshold_sq.sqrt().acos().to_degrees(),
        slow_dihedral.acos().to_degrees()
    );

    let threshold_sq = 54.0_f32.to_radians().cos().powi(2);
    let result =
        IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_same_winding(a, b, c, d, threshold_sq);
    assert_eq!(
        result,
        None,
        "failed cos at {}° slow_dihedral:{}°",
        threshold_sq.sqrt().acos().to_degrees(),
        slow_dihedral.acos().to_degrees()
    );
}

#[test]
fn test_arrowhead_quad_3() {
    let a: Vec3A = [-0.25, -1.0, -0.75].into();
    let b: Vec3A = [-0.5, -1.0, -1.0].into();
    let c: Vec3A = [0.5, -1.0, -1.0].into();
    let d: Vec3A = [-0.125, -1.0, -0.375].into();

    let threshold = 0.9;
    let result = IsotropicRemeshAlgo::<f32, Vec3>::quad_is_convex_coplanar(a, b, c, d, threshold);
    assert!(!result);
}

#[test]
fn test_some_triangles() {
    let threshold_sq = 0.985_f32.powi(2);

    let v = [
        [-1.3919225, -0.38369077, 1.2062045],
        [1.3399366, -0.38369077, 1.2062045],
        [1.3399366, 0.29927406, 0.4109751],
        [-1.3919225, -1.1789197, 1.3684515],
    ];
    let a = v[0].into();
    let b = v[1].into();
    let c = v[2].into();
    let d = v[3].into();

    let result =
        IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_same_winding(a, b, c, d, threshold_sq);
    assert_eq!(result, None);

    let v = [
        [-1.3919225, -0.38369077, 1.2062045],
        [1.3399366, 0.29927406, 0.4109751],
        [-0.025992602, 0.79981124, -0.2719897],
        [-1.3919225, -1.1789197, 1.3684515],
    ];
    let a = v[0].into();
    let b = v[1].into();
    let c = v[2].into();
    let d = v[3].into();

    let result =
        IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_same_winding(a, b, c, d, threshold_sq);
    assert_eq!(result, None);

    let v = [
        [-1.3919225, -0.38369077, 1.2062045],
        [-0.025992602, 0.79981124, -0.2719897],
        [-1.3919225, 0.79981124, -0.2719897],
        [-1.3919225, -1.1789197, 1.3684515],
    ];
    let a = v[0].into();
    let b = v[1].into();
    let c = v[2].into();
    let d = v[3].into();

    let result =
        IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_same_winding(a, b, c, d, threshold_sq);
    assert_eq!(result, None);

    let v = [
        [-1.3919225, -0.38369077, 1.2062045],
        [1.3399366, -1.1789197, 1.3684515],
        [1.3399366, -0.38369077, 1.2062045],
        [-1.3919225, -1.1789197, 1.3684515],
    ];
    let a = v[0].into();
    let b = v[1].into();
    let c = v[2].into();
    let d = v[3].into();

    let result =
        IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_same_winding(a, b, c, d, threshold_sq);
    assert_eq!(result, Some(true));
}

#[test]
fn test_normal_similarity_for_collapse_rotated() {
    let threshold_sq = 0.9_f32.powi(2);

    // Test rotation around Z axis at various angles
    for angle_deg in [0.0_f32, 15.0, 30.0, 45.0, 60.0, 75.0, 90.0] {
        let angle = angle_deg.to_radians();
        let rotation = Mat3::from_rotation_z(angle);

        // Original coplanar triangles in XY plane
        let a = rotation * Vec3A::new(0.0, 0.0, 0.0);
        let b = rotation * Vec3A::new(0.0, 1.0, 0.0);
        let c = rotation * Vec3A::new(1.0, 1.0, 0.0);
        let d = rotation * Vec3A::new(1.0, 0.0, 0.0);

        let result = IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_same_winding(
            a,
            b,
            c,
            d,
            threshold_sq,
        );
        assert_eq!(
            result,
            Some(true),
            "Failed at {}° rotation around Z",
            angle_deg
        );
    }

    // Also test rotation around X axis
    for angle_deg in [0.0_f32, 15.0, 30.0, 45.0, 60.0, 75.0, 90.0] {
        let angle = angle_deg.to_radians();
        let rotation = Mat3::from_rotation_x(angle);

        let a = rotation * Vec3A::new(0.0, 0.0, 0.0);
        let b = rotation * Vec3A::new(0.0, 1.0, 0.0);
        let c = rotation * Vec3A::new(1.0, 1.0, 0.0);
        let d = rotation * Vec3A::new(1.0, 0.0, 0.0);

        let result = IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_same_winding(
            a,
            b,
            c,
            d,
            threshold_sq,
        );
        assert_eq!(
            result,
            Some(true),
            "Failed at {} degrees rotation around X",
            angle_deg
        );
    }

    // And rotation around Y axis
    for angle_deg in [0.0_f32, 15.0, 30.0, 45.0, 60.0, 75.0, 90.0] {
        let angle = angle_deg.to_radians();
        let rotation = Mat3::from_rotation_y(angle);

        let a = rotation * Vec3A::new(0.0, 0.0, 0.0);
        let b = rotation * Vec3A::new(0.0, 1.0, 0.0);
        let c = rotation * Vec3A::new(1.0, 1.0, 0.0);
        let d = rotation * Vec3A::new(1.0, 0.0, 0.0);

        let result = IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_same_winding(
            a,
            b,
            c,
            d,
            threshold_sq,
        );
        assert_eq!(
            result,
            Some(true),
            "Failed at {} degrees rotation around Y",
            angle_deg
        );
    }
}

#[test]
fn test_triangles_not_inverted_1() {
    // Simple coplanar triangles sharing edge BC
    // Triangle ABC in xy-plane with normal pointing up (+z)
    let va = Vec3A::new(0.0, 0.0, 0.0); // A at origin
    let vb = Vec3A::new(1.0, 0.0, 0.0); // B on x-axis
    let vc = Vec3A::new(0.0, 1.0, 0.0); // C on y-axis

    // Triangle DBC with D creating another triangle in same plane
    let vd = Vec3A::new(1.0, 1.0, 0.0); // D completes a square

    let cos_threshold_sq = -0.9 * -0.9;

    let result = IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_opposite_winding(
        va,
        vb,
        vc,
        vd,
        cos_threshold_sq,
    );

    // Both triangles have normals pointing in +z direction
    // They should be compatible (not inverted)
    assert_eq!(result, Some(true));
}

#[test]
fn test_triangles_same_orientation() {
    // Two triangles in same plane with normals pointing same direction
    let va = Vec3A::new(0.0, 0.0, 0.0);
    let vb = Vec3A::new(1.0, 0.0, 0.0);
    let vc = Vec3A::new(0.0, 1.0, 0.0);
    let vd = Vec3A::new(1.0, 1.0, 0.0);

    let result = IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_opposite_winding(
        va,
        vb,
        vc,
        vd,
        -0.9 * -0.9,
    );

    assert_eq!(result, Some(true));
}

#[test]
fn test_degenerate_triangle_abc() {
    // Triangle ABC is degenerate - A, B, C are collinear
    let va = Vec3A::new(0.0, 0.0, 0.0);
    let vb = Vec3A::new(1.0, 0.0, 0.0);
    let vc = Vec3A::new(2.0, 0.0, 0.0); // Collinear with A and B

    let vd = Vec3A::new(1.0, 1.0, 0.0);

    let result = IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_opposite_winding(
        va,
        vb,
        vc,
        vd,
        -0.9 * -0.9,
    );

    assert_eq!(result, Some(false)); // ABC is degenerate
}

#[test]
fn test_nearly_inverted_at_threshold() {
    // Triangle ABC pointing up
    let va = Vec3A::new(0.0, 0.0, 0.0);
    let vb = Vec3A::new(1.0, 0.0, 0.0);
    let vc = Vec3A::new(0.0, 1.0, 0.0);

    // D positioned to create ~154° angle (cos ≈ -0.9)
    // At 154°, the triangles are nearly opposite but just acceptable
    let vd = Vec3A::new(1.0, 1.0, -0.4);

    let result = IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_opposite_winding(
        va,
        vb,
        vc,
        vd,
        -0.9 * -0.9,
    );

    // Should be close to threshold - might be Some(true) or None depending on exact value
    assert!(result.is_some() || result.is_none());
}

#[test]
fn test_obtuse_angle_acceptable() {
    // Triangle ABC pointing up
    let va = Vec3A::new(0.0, 0.0, 0.0);
    let vb = Vec3A::new(1.0, 0.0, 0.0);
    let vc = Vec3A::new(0.0, 1.0, 0.0);

    // D creates an obtuse angle (~120°) but not inverted
    let vd = Vec3A::new(1.0, 1.0, -0.1);

    let result = IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_opposite_winding(
        va,
        vb,
        vc,
        vd,
        -0.9 * -0.9,
    );

    // cos(120°) = -0.5, which is > -0.9, so should pass
    assert_eq!(result, Some(true));
}

#[test]
fn test_strict_threshold() {
    // Triangle ABC pointing up
    let va = Vec3A::new(0.0, 0.0, 0.0);
    let vb = Vec3A::new(1.0, 0.0, 0.0);
    let vc = Vec3A::new(0.0, 1.0, 0.0);

    // D creates ~135° angle
    let vd = Vec3A::new(1.0, 1.0, -0.2);

    // With stricter threshold (-0.5 = 120°), this should pass
    let result = IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_opposite_winding(
        va,
        vb,
        vc,
        vd,
        -0.5 * -0.5,
    );
    assert_eq!(result, Some(true));

    // With very permissive threshold (-0.95 = 162°), this should also pass
    let result = IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_opposite_winding(
        va,
        vb,
        vc,
        vd,
        -0.95 * -0.95,
    );
    assert_eq!(result, Some(true));
}

#[test]
fn test_normal_similarity_for_collapse() {
    // Triangle ABC pointing up
    let va: Vec3A = [-0.5, -0.5, 0.5].into();
    let vb: Vec3A = [0.5, -0.5, -0.5].into();
    let vc: Vec3A = [0.5, -0.5, 0.5].into();
    let vd: Vec3A = [-0.5, -0.5, -0.5].into();

    let threshold = std::f32::consts::FRAC_1_SQRT_2;

    let dot_slow =
        IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_same_winding_slow(va, vb, vc, vd);
    let dot_fast = IsotropicRemeshAlgo::<f32, Vec3>::normal_similarity_same_winding(
        va,
        vb,
        vc,
        vd,
        threshold * threshold,
    );
    println!("dot_slow = {}", dot_slow);
    println!("dot_fast = {:?}", dot_fast);

    assert!(dot_slow > threshold);
    assert_eq!(dot_fast, Some(true));
}

#[test]
fn test_corner_pool() -> Result<(), RemeshError> {
    let vertices = [
        [-1.0, -1.0, -1.0],
        [-1.0, -1.0, 1.0],
        [-1.0, 1.0, -1.0],
        [-1.0, 1.0, 1.0],
        [1.0, -1.0, -1.0],
        [1.0, -1.0, 1.0],
        [1.0, 1.0, -1.0],
        [1.0, 1.0, 1.0],
    ];
    let indices = [
        1, 2, 0, 3, 6, 2, 7, 4, 6, 5, 0, 4, 6, 0, 2, 3, 5, 7, 1, 3, 2, 3, 7, 6, 7, 5, 4, 5, 1, 0,
        6, 4, 0, 3, 1, 5,
    ];

    let remesher = IsotropicRemesh::<f32, _, true>::new(&vertices, &indices)?
        .with_target_edge_length(0.5_f32)?
        .with_default_collapse_multiplier()?
        .with_default_split_multiplier()?
        .with_default_flip_edges()?
        .with_smooth_weight(0.01)?
        .with_print_stats(10)?
        .with_flip_edges(FlipStrategy::Valence)?;

    let mut algo = remesher.new_algo()?;
    algo.run_internal(10)?;

    assert_eq!(
        algo.vertex_pool.active_count(),
        142,
        "result_vertices.len() mismatch"
    );
    assert_eq!(
        algo.corner_table.triangle_pool.active_count(),
        280,
        "triangles.len() mismatch"
    );

    algo.with_target_edge_length(2.2)?;
    algo.run_internal(20)?;
    let (result_vertices, result_indices) = algo.compress_mesh()?;

    assert_eq!(result_vertices.len(), 8, "result_vertices.len() mismatch");
    assert_eq!(result_indices.len(), 36, "result_indices.len() mismatch");
    Ok(())
}

#[test]
fn test_debug_impl() {
    // fill up code coverage.
    let _ = format!("{:?}", VertexIndex(0));
    let _ = format!("{:?}", CornerIndex(0));
    let _ = format!("{:?}", TriangleIndex(0));
    let _ = format!("{:?}", VertexIndex::INVALID);
    let _ = format!("{:?}", CornerIndex::INVALID);
    let _ = format!("{:?}", TriangleIndex::INVALID);
}