// example_flatten.cpp
//
// conformallab++ — Real conformal flattening example
//
// This example demonstrates the PRIMARY USE CASE of the library:
// conformally flatten a 3-D surface mesh to the plane with minimal
// angle distortion.  Every interior vertex is assigned a target cone
// angle of 2π (a regular flat vertex); the solver finds the unique
// conformal factor u_v that realises this target.
//
// Contrast with the other examples (example_euclidean, example_layout)
// which use the "natural theta" testing trick that produces x* = 0 —
// a valid solver test but NOT a conformal flattening.
//
// ┌─────────────────────────────────────────────────────────────────────┐
// │  Pipeline for conformal flattening of an OPEN mesh (disk topology) │
// │                                                                     │
// │  1. Load mesh                                                       │
// │  2. Setup maps + compute λ° from geometry                          │
// │  3. Pin all boundary vertices (they define the boundary of the UV) │
// │     Set Θ_v = 2π for all interior vertices (flat target)           │
// │     ← NO Gauss–Bonnet check needed for open meshes                 │
// │  4. Solve Newton                                                    │
// │  5. Compute planar layout — the conformal UV parameterisation       │
// │  6. Save UV-mapped OFF + report distortion                         │
// └─────────────────────────────────────────────────────────────────────┘
//
// Build (requires -DWITH_CGAL=ON):
//   cmake -S code -B build -DWITH_CGAL=ON
//   cmake --build build --target example_flatten
//
// Run:
//   # With a real 3-D surface mesh (open, disk topology):
//   ./build/examples/example_flatten code/data/obj/cathead.obj flat.off
//
//   # Without arguments: uses a built-in synthetic open mesh (6 vertices)
//   ./build/examples/example_flatten

#include "conformal_mesh.hpp"
#include "mesh_builder.hpp"
#include "mesh_io.hpp"
#include "euclidean_functional.hpp"
#include "gauss_bonnet.hpp"
#include "newton_solver.hpp"
#include "layout.hpp"
#include "constants.hpp"
#include <CGAL/boost/graph/iterator.h>
#include <iostream>
#include <iomanip>
#include <string>
#include <vector>
#include <cmath>
#include <algorithm>

using namespace conformallab;

int main(int argc, char* argv[])
{
    // ── Step 1: load or synthesise a mesh ────────────────────────────────────
    ConformalMesh mesh;
    std::string input_path  = (argc > 1) ? argv[1] : "";
    std::string output_path = (argc > 2) ? argv[2] : "/tmp/conformallab_flatten_out.off";

    if (input_path.empty()) {
        // Fallback: load cathead.obj from the standard data location if
        // it exists alongside the executable; otherwise use a synthetic mesh.
        // For a meaningful non-trivial flattening, supply a real 3-D mesh:
        //   ./example_flatten code/data/obj/cathead.obj
        std::cout << "[example_flatten] No input given.  Using make_quad_strip() "
                     "(flat synthetic mesh — u_v will be near-zero).\n"
                  << "  For a non-trivial flattening, provide a 3-D mesh:\n"
                  << "    ./example_flatten code/data/obj/cathead.obj\n\n";
        mesh = make_quad_strip();
    } else {
        std::cout << "[example_flatten] Loading mesh: " << input_path << "\n";
        try { mesh = load_mesh(input_path); }
        catch (const std::exception& e) {
            std::cerr << "Error loading mesh: " << e.what() << "\n";
            return 1;
        }
    }

    const int V = static_cast<int>(mesh.number_of_vertices());
    const int F = static_cast<int>(mesh.number_of_faces());
    std::cout << "[example_flatten] Mesh: " << V << " vertices, " << F << " faces\n";

    // ── Step 2: setup maps + compute initial edge lengths ─────────────────────
    auto maps = setup_euclidean_maps(mesh);
    compute_euclidean_lambda0_from_mesh(mesh, maps);

    // ── Step 3: DOF assignment — pin boundary, free interior ──────────────────
    //
    // For an OPEN mesh (disk topology):
    //   • Boundary vertices are pinned (v_idx = -1): they define the
    //     boundary of the UV domain and are not optimised.
    //   • Interior vertices get a DOF (v_idx ≥ 0) and target Θ_v = 2π.
    //     This means "make every interior point look like a flat plane vertex".
    //
    // For a CLOSED mesh (e.g. a sphere or torus), use the Euclidean pipeline
    // on a cut-open mesh (see example_layout.cpp + cut_graph.hpp), or use the
    // spherical / hyper-ideal functional instead.
    //
    // ⚠  This is the KEY difference from example_euclidean.cpp:
    //    There, "natural theta" sets Θ_v = actual angle sum → x* = 0 (trivial).
    //    Here, Θ_v = 2π → the solver finds the REAL conformal map.
    int idx = 0;
    int n_boundary = 0, n_interior = 0;

    for (auto v : mesh.vertices()) {
        // CGAL: a vertex is on the boundary iff it has an incident border halfedge.
        bool is_bnd = false;
        for (auto h : CGAL::halfedges_around_target(v, mesh))
            if (mesh.is_border(h)) { is_bnd = true; break; }

        if (is_bnd) {
            maps.v_idx[v] = -1;      // pinned — boundary defines the UV border
            ++n_boundary;
        } else {
            maps.theta_v[v] = TWO_PI;  // flat interior target — the actual goal
            maps.v_idx[v]   = idx++;
            ++n_interior;
        }
    }

    std::cout << "[example_flatten] Boundary vertices (pinned): " << n_boundary
              << "  Interior (free DOFs): " << n_interior << "\n";

    if (n_interior == 0) {
        std::cerr << "[example_flatten] No interior vertices — mesh has no free DOFs.\n"
                  << "  Use a mesh with interior vertices (e.g. cathead.obj).\n";
        return 1;
    }

    // For open meshes the Gauss–Bonnet identity holds in a different form and
    // does NOT need to be checked before calling newton_euclidean.  The solver
    // converges as long as at least one interior vertex exists.
    // (For CLOSED meshes: call enforce_gauss_bonnet(mesh, maps) here.)

    // ── Step 4: Newton ────────────────────────────────────────────────────────
    std::vector<double> x0(static_cast<std::size_t>(idx), 0.0);

    std::cout << "[example_flatten] Running Newton (tol = 1e-9)…\n";
    auto res = newton_euclidean(mesh, x0, maps, /*tol=*/1e-9);

    if (res.converged)
        std::cout << "[example_flatten] Converged in " << res.iterations
                  << " iterations  ||G||_inf = " << std::scientific
                  << std::setprecision(2) << res.grad_inf_norm << "\n";
    else
        std::cout << "[example_flatten] WARNING: did not converge after "
                  << res.iterations << " iterations  ||G||_inf = "
                  << res.grad_inf_norm << "\n";

    // ── Step 5: report conformal factors ──────────────────────────────────────
    //
    // u_v is the log-scale factor: the area element at vertex v is scaled by
    // exp(2·u_v).  For a non-trivial mesh the values are non-zero.
    double u_min = 0.0, u_max = 0.0;
    for (auto v : mesh.vertices()) {
        int iv = maps.v_idx[v];
        if (iv >= 0) {
            double u = res.x[static_cast<std::size_t>(iv)];
            u_min = std::min(u_min, u);
            u_max = std::max(u_max, u);
        }
    }
    std::cout << "[example_flatten] Conformal factors u_v: "
              << "min = " << std::fixed << std::setprecision(4) << u_min
              << "  max = " << u_max << "\n";

    if (std::abs(u_max - u_min) < 1e-6)
        std::cout << "[example_flatten] Note: u_v ≈ 0 everywhere — the mesh is "
                     "already flat in the\n"
                     "  Euclidean sense (e.g. a planar mesh).  Use a 3-D surface "
                     "for non-trivial output.\n";
    else
        std::cout << "[example_flatten] Non-trivial u_v range = "
                  << (u_max - u_min) << " — real conformal deformation computed.\n";

    // ── Step 6: compute and save UV layout ────────────────────────────────────
    auto layout = euclidean_layout(mesh, res.x, maps);

    if (layout.success) {
        // Find the bounding box of the UV coords
        double xmin = 1e30, xmax = -1e30, ymin = 1e30, ymax = -1e30;
        for (auto v : mesh.vertices()) {
            if (static_cast<std::size_t>(v.idx()) < layout.uv.size()) {
                const auto& p = layout.uv[static_cast<std::size_t>(v.idx())];
                xmin = std::min(xmin, p.x()); xmax = std::max(xmax, p.x());
                ymin = std::min(ymin, p.y()); ymax = std::max(ymax, p.y());
            }
        }
        std::cout << "[example_flatten] UV bounding box: ["
                  << std::fixed << std::setprecision(3)
                  << xmin << ", " << xmax << "] × ["
                  << ymin << ", " << ymax << "]\n";

        save_layout_off(output_path, mesh, layout);
        std::cout << "[example_flatten] UV layout saved → " << output_path << "\n"
                  << "  Open in MeshLab or Blender to inspect the UV parameterisation.\n";
    } else {
        std::cerr << "[example_flatten] Layout failed.\n";
    }

    return (res.converged && layout.success) ? 0 : 1;
}