main.cpp

#include "step1.cpp"
#include "transforms.cpp"
#include <algorithm>
#include <cmath>
#include <iostream>
#include <numeric>
#include <iomanip>
#include <string>
#include <filesystem>
#include <iostream>
#include "mmg/mmg2d/libmmg2d.h"

void init_directories() {
    namespace fs = std::filesystem;
    fs::create_directory("../temp");
    fs::create_directory("../tests_in");
    fs::create_directories("../tests_out/side_by_side");
}


std::string tests[] = {
    "shift", 
    "full_isometry",
    "full_repulsion",
    "simple_rotation", 
    "full_rotation",
    "vertical_to_horizontal",
};


void vertical_to_horizontal() {
    PairedTriangles m {};
    auto attributes = m.import_raw("../tests_base/vertical_to_horizontal.geogram");
    m.resize_attrs();
    PointAttribute<bool> selection("selection", attributes, m);
    m.connect();
    simple_rotation(m, selection, {.5, .5, 0}, -M_PI_2);
    m.save("../tests_in/vertical_to_horizontal");
}

void full_isometry() {
    PairedTriangles m {};
    auto attributes = m.import_raw("../tests_base/full_rotation.geogram");
    m.resize_attrs();
    PointAttribute<bool> selection("selection", attributes, m);
    m.connect();
    simple_isometry(m, selection, {.5, .4, 0}, 1.7);
    m.save("../tests_in/full_isometry");
}

void full_repulsion() {
    PairedTriangles m {};
    auto attributes = m.import_raw("../tests_base/full_rotation.geogram");
    m.resize_attrs();
    PointAttribute<bool> selection("selection", attributes, m);
    m.connect();
    local_repulsion(m, selection, {.4, .5, 0}, 0.4, 0.12);
    m.save("../tests_in/full_repulsion");
}

void simple_rotation() {
    PairedTriangles m {};
    auto attributes = m.import_raw("../tests_base/full_rotation.geogram");
    m.resize_attrs();
    PointAttribute<bool> selection("selection", attributes, m);
    m.connect();
    simple_rotation(m, selection, {.5, .5, 0}, -M_PI);
    m.save("../tests_in/simple_rotation");
}

void full_rotation() {
    PairedTriangles m {};
    auto attributes = m.import_raw("../tests_base/full_rotation.geogram");
    m.resize_attrs();
    PointAttribute<bool> selection("selection", attributes, m);
    m.connect();
    simple_shift(m, selection, 0., 0.2);
    simple_rotation(m, selection, {.5, .5, 0}, -M_PI);
    m.save("../tests_in/full_rotation");
}

void full_rotation_test() {
    PairedTriangles m {};
    auto attributes = m.import_raw("../tests_base/full_rotation_test.geogram");
    m.resize_attrs();
    PointAttribute<bool> selection("selection", attributes, m);
    m.connect();
    simple_shift(m, selection, 0.2, 0.);
    simple_rotation(m, selection, {.5, .5, 0}, -2.*M_PI);
    m.save("../tests_in/full_rotation_test");
}

void shift() {
    PairedTriangles m {};
    auto attributes = m.import_raw("../tests_base/shift.geogram");
    m.resize_attrs();
    PointAttribute<bool> selection("selection", attributes, m);
    m.connect();
    simple_shift(m, selection, 0.2, 0.);
    m.save("../tests_in/shift");
}

void init_transforms() {
    vertical_to_horizontal();
    full_isometry();
    full_repulsion();
    simple_rotation();
    full_rotation();
    shift();
}


void run_tests() {
    srand(time(NULL));
    for (std::string fname: tests) {
        std::cout << "---------------------------------\n";
        std::cout << "START " + fname + "\n\n\n";
        // INIT

        PairedTriangles dummy;
        dummy.import("../tests_in/" + fname);
        dummy.connect();
        dummy.resize_attrs();
        dummy.dyn_facets = dummy.nfacets();
        std::cout << dummy.dyn_facets << "\n";
        for (auto f: dummy.iter_facets()) {
            dummy.energy[f] = 0;
            dummy.update(f);
        }
        std::cout << "before:" << dummy.mean_energy() << "\n";
        const auto temp_fname = "../temp/" + fname + "_remeshed";
        const auto temp_dir = temp_fname + ".geogram";
        dummy.save(temp_fname);
        dummy.save_animation_path("../tests_out/side_by_side/" + fname + "_before.geogram");

        int N = 5;
        
        for (int i: range(N)) {

            std::cout << fname << ": " << i << "\n";

            PairedTriangles m {};
            m.import(temp_fname);
            m.connect();
            m.resize_attrs();
            m.dyn_facets = m.nfacets();
            for (auto f: m.iter_facets()) {
                m.energy[f] = 0;
                m.update(f);
            }
            run(m, fname);
            m.save(temp_fname);
            m.save("../tests_out/"+fname);

            if (i == N-1) break;

            MMG5_pMesh mesh = NULL;
            MMG5_pSol met = NULL;
            m.to_mmg(mesh, met);
            int ier = MMG2D_mmg2dlib(mesh, met);
            if (ier == MMG5_STRONGFAILURE) {
                throw std::runtime_error("MMG severe error: remeshing failed.");
            } else if (ier == MMG5_LOWFAILURE) {
                std::cerr << "MMG warning: remeshing completed with low-priority errors.\n";
            }
            std::string bucket_dir = "../temp/bucket.mesh";
            MMG2D_saveMesh(mesh, bucket_dir.c_str());
            PairedTriangles remeshed {};
            remeshed.import_raw(bucket_dir);
            remeshed.resize_attrs();
            remeshed.connect();
            double eps = 1e-14;
            
            for (auto v: remeshed.iter_vertices()) {
                remeshed.I[v].z = 0;
                remeshed.J[v].z = 0;
                for (auto f: m.iter_facets()) { // original facets
                    vec3 b = map(m.I, f).bary_coords(remeshed.I[v]);
                    if (b.x >= -eps && b.y >= -eps && b.z >= -eps) {
                        remeshed.J[v] = ((map(m.J, f).as_matrix()).transpose() * b);
                    }
                }
            }
            untangle_symmetric(remeshed);
        }
        
        PairedTriangles m {};
        m.import("../tests_out/" + fname);
        m.save_animation_path("../tests_out/side_by_side/" + fname + "_after.geogram");
    }
}

struct Stats { // CHADGPT
    double mean, sigma, min, max;
};

auto compute_stats = [](const auto& v) -> Stats {
    const double mean =
        std::accumulate(v.begin(), v.end(), 0.0) / v.size();

    double variance = 0.0;
    for (double x : v)
        variance += (x - mean) * (x - mean);
    variance /= v.size();

    const auto [min_it, max_it] =
        std::minmax_element(v.begin(), v.end());

    return {
        mean,
        std::sqrt(variance),
        *min_it,
        *max_it
    };
};

void dump_tests() {
    for (auto fname: tests) {
        PairedTriangles before, after, mmg;
        before.import("../tests_in/" + fname);
        after.import("../tests_out/" + fname);

        before.resize_attrs();
        after.resize_attrs();

        std::cout << "---------------" + fname + "-------------------\n\n";

        for (auto f: before.iter_facets()) {
            before.update(f);
        }
        for (auto f: after.iter_facets()) {
            after.update(f);
        }

        {
            auto b = compute_stats(before.energy.ptr->data);
            auto a = compute_stats(after.energy.ptr->data);

            std::cout << "energy: \n\n";
            std::cout
                << std::setw(12) << ""      << std::setw(15) << "before" << std::setw(15) << "after" << '\n'
                << std::setw(12) << "mean"  << std::setw(15) << b.mean   << std::setw(15) << a.mean  << '\n'
                << std::setw(12) << "sigma" << std::setw(15) << b.sigma  << std::setw(15) << a.sigma << '\n'
                << std::setw(12) << "min"   << std::setw(15) << b.min    << std::setw(15) << a.min   << '\n'
                << std::setw(12) << "max"   << std::setw(15) << b.max    << std::setw(15) << a.max   << '\n';
        }

        {
            auto b = compute_stats(before.angles());
            auto a = compute_stats(after.angles());

            std::cout << "\n\nangles (degrees): \n\n";
            std::cout
                << std::setw(12) << ""      << std::setw(15) << "before" << std::setw(15) << "after" << '\n'
                << std::setw(12) << "mean"  << std::setw(15) << b.mean   << std::setw(15) << a.mean  << '\n'
                << std::setw(12) << "sigma" << std::setw(15) << b.sigma  << std::setw(15) << a.sigma << '\n'
                << std::setw(12) << "min"   << std::setw(15) << b.min    << std::setw(15) << a.min   << '\n'
                << std::setw(12) << "max"   << std::setw(15) << b.max    << std::setw(15) << a.max   << '\n';
        }

        {
            auto b = compute_stats(before.aspect_ratios());
            auto a = compute_stats(after.aspect_ratios());

            std::cout << "\n\naspect ratios: \n\n";
            std::cout
                << std::setw(12) << ""      << std::setw(15) << "before" << std::setw(15) << "after" << '\n'
                << std::setw(12) << "mean"  << std::setw(15) << b.mean   << std::setw(15) << a.mean  << '\n'
                << std::setw(12) << "sigma" << std::setw(15) << b.sigma  << std::setw(15) << a.sigma << '\n'
                << std::setw(12) << "min"   << std::setw(15) << b.min    << std::setw(15) << a.min   << '\n'
                << std::setw(12) << "max"   << std::setw(15) << b.max    << std::setw(15) << a.max   << '\n';
        }

        {
            auto b = compute_stats(before.distortion());
            auto a = compute_stats(after.distortion());

            std::cout << "\n\ndistortion: \n\n";
            std::cout
                << std::setw(12) << ""      << std::setw(15) << "before" << std::setw(15) << "after" << '\n'
                << std::setw(12) << "mean"  << std::setw(15) << b.mean   << std::setw(15) << a.mean  << '\n'
                << std::setw(12) << "sigma" << std::setw(15) << b.sigma  << std::setw(15) << a.sigma << '\n'
                << std::setw(12) << "min"   << std::setw(15) << b.min    << std::setw(15) << a.min   << '\n'
                << std::setw(12) << "max"   << std::setw(15) << b.max    << std::setw(15) << a.max   << '\n';
        }
    }
}

int main() {
    std::cout << "Initializing directories" << "\n";
    init_directories();
    std::cout << "Generating inputs" << "\n";
    init_transforms();
    std::cout << "Building outputs" << "\n";
    run_tests();
    std::cout << "Computing results" << "\n";
    dump_tests();
    std::cout << "Done" << "\n";
    return 0;
}