diff --git a/app/qslim/app/main.cpp b/app/qslim/app/main.cpp index b6eb90a5c7..9f370186bc 100644 --- a/app/qslim/app/main.cpp +++ b/app/qslim/app/main.cpp @@ -99,7 +99,6 @@ int main(int argc, char** argv) const double diag = (box_minmax.first - box_minmax.second).norm(); const double envelope_size = env_rel * diag; - igl::Timer timer; timer.start(); QSLIM m(verts, num_thread); diff --git a/app/remeshing/app/main.cpp b/app/remeshing/app/main.cpp index 8f6d3a4f2a..f4bec2c9a8 100644 --- a/app/remeshing/app/main.cpp +++ b/app/remeshing/app/main.cpp @@ -9,7 +9,6 @@ #include #include - #include #include #include @@ -23,26 +22,41 @@ using namespace std::chrono; #include "remeshing_spec.hpp" -void run_remeshing(std::string input, double len, std::string output, UniformRemeshing& m, int itrs) +void run_remeshing( + std::string input, + std::string output, + UniformRemeshing& m, + int itrs, + bool debug_output = false) { + // wmtk::logger().info( + // "Before_vertices#: {} \n Before_tris#: {}", + // m.vert_capacity(), + // m.tri_capacity()); + auto start = high_resolution_clock::now(); - wmtk::logger().info("target len: {}", len); - m.uniform_remeshing(len, itrs); + m.uniform_remeshing(itrs, debug_output); // m.consolidate_mesh(); auto stop = high_resolution_clock::now(); auto duration = duration_cast(stop - start); m.consolidate_mesh(); - m.write_triangle_mesh(output); + + std::filesystem::path outp(output); + std::string base = outp.replace_extension("").string(); + m.write_vtu(output); + + // m.write_triangle_mesh(output); + // m.write_feature_vertices_obj(output + ".feature_vertices.obj"); auto properties = m.average_len_valen(); - wmtk::logger().info("runtime in ms {}", duration.count()); - wmtk::logger().info("current_memory {}", getCurrentRSS() / (1024. * 1024)); - wmtk::logger().info("peak_memory {}", getPeakRSS() / (1024. * 1024)); - wmtk::logger().info("after remesh properties: {}", properties); - wmtk::logger().info( - "After_vertices#: {} \n\t After_tris#: {}", - m.vert_capacity(), - m.tri_capacity()); + // wmtk::logger().info("runtime in ms {}", duration.count()); + // wmtk::logger().info("current_memory {}", getCurrentRSS() / (1024. * 1024)); + // wmtk::logger().info("peak_memory {}", getPeakRSS() / (1024. * 1024)); + // wmtk::logger().info("after remesh properties: {}", properties); + // wmtk::logger().info( + // "After_vertices#: {} \n\t After_tris#: {}", + // m.vert_capacity(), + // m.tri_capacity()); } int main(int argc, char** argv) @@ -65,8 +79,7 @@ int main(int argc, char** argv) std::ifstream ifs(json_input_file); json_params = nlohmann::json::parse(ifs); } catch (const std::exception& e) { - logger().error("Could not load or parse JSON input file"); - logger().error(e.what()); + log_and_throw_error("Could not load or parse JSON input file: {}", e.what()); } // verify input and inject defaults @@ -79,6 +92,79 @@ int main(int argc, char** argv) json_params = spec_engine.inject_defaults(json_params, remeshing_spec); } + + std::vector> fixed_vertices; + std::vector, int>> feature_edge_list; + + std::vector face_patch_ids; + + try { + if (json_params.contains("patches")) { + std::string patches_path = json_params["patches"]; + logger().info("Loading patches info file: {}", patches_path); + + + nlohmann::json patches_json; + std::ifstream cifs(patches_path); + if (!cifs.is_open()) { + logger().error("Cannot open patches file {}", patches_path); + } else { + cifs >> patches_json; + } + + const auto& corners = patches_json["corner_vids"]; + fixed_vertices.reserve(corners.size()); + + for (const auto& [key, val] : corners.items()) { + const int corner_ids = std::stoi(key) + 1; // add 1 so 0 means "not frozen" + const size_t vid = val.get(); + + fixed_vertices.emplace_back(vid, corner_ids); + } + logger().info("Loaded {} fixed (corner) vertices.", fixed_vertices.size()); + + const auto& feature_edges = patches_json["edge2seg"]; + feature_edge_list.reserve(feature_edges.size()); + + for (const auto& [key, val] : feature_edges.items()) { + const auto comma_pos = key.find(','); + if (comma_pos == std::string::npos) { + logger().error("Invalid feature edge key (no comma): {}", key); + continue; + } + + const size_t v0 = std::stoul(key.substr(0, comma_pos)); + const size_t v1 = std::stoul(key.substr(comma_pos + 1)); + + const int seg_id_plus1 = val.get() + 1; // <-- THIS is the id (+1) + + feature_edge_list.push_back({{{v0, v1}}, seg_id_plus1}); + } + + logger().info("Loaded {} feature edges.", feature_edge_list.size()); + + const auto& pids = patches_json["fid2patch"]; + face_patch_ids.resize(pids.size()); + for (const auto& [key, pid] : pids.items()) { + const size_t fid = std::stoul(key); + face_patch_ids[fid] = pid; + } + } + } catch (const std::exception& e) { + log_and_throw_error("Error reading corner vertices: {}", e.what()); + } + + std::vector feature_vertices; + feature_vertices.reserve(feature_edge_list.size() * 2); + + for (const auto& e : feature_edge_list) { + const auto& ab = e.first; + feature_vertices.push_back(ab[0]); + feature_vertices.push_back(ab[1]); + } + + wmtk::vector_unique(feature_vertices); + // logger settings { std::string log_file_name = json_params["log_file"]; @@ -90,37 +176,87 @@ int main(int argc, char** argv) const std::string input_path = json_params["input"]; const std::string output = json_params["output"]; - const double env_rel = json_params["eps_rel"]; + double eps = json_params["eps"]; + const double eps_rel = json_params["eps_rel"]; const double length_rel = json_params["length_rel"]; + const double length_factor = json_params["length_factor"]; const int num_threads = json_params["num_threads"]; const int itrs = json_params["max_iterations"]; const bool sample_envelope = json_params["use_sample_envelope"]; const bool freeze_boundary = json_params["freeze_boundary"]; double length_abs = json_params["length_abs"]; + const bool debug_output = json_params["DEBUG_output"]; wmtk::logger().info("remeshing on {}", input_path); wmtk::logger().info("freeze bnd {}", freeze_boundary); + // create report file + auto report_file = output + "_report.json"; + + // Avoid confusion: if the report already exists, delete it and write a fresh one. + if (std::filesystem::exists(report_file)) { + std::filesystem::remove(report_file); + } + + nlohmann::json report = nlohmann::json::object(); + report["before"] = nlohmann::json::object(); + report["after"] = nlohmann::json::object(); std::vector verts; std::vector> tris; std::pair box_minmax; - double remove_duplicate_eps = 1e-5; - std::vector modified_nonmanifold_v; - wmtk::stl_to_manifold_wmtk_input( - input_path, - remove_duplicate_eps, - box_minmax, - verts, - tris, - modified_nonmanifold_v); + Eigen::MatrixXd inV; + Eigen::MatrixXi inF; + igl::read_triangle_mesh(input_path, inV, inF); + verts.resize(inV.rows()); + tris.resize(inF.rows()); + wmtk::eigen_to_wmtk_input(verts, tris, inV, inF); + + { + // basic mesh size stats + report["before"]["nV"] = static_cast(inV.rows()); + report["before"]["nF"] = static_cast(inF.rows()); + + // min/max internal angle + { + Eigen::VectorXd angles; + igl::internal_angles(inV, inF, angles); + auto min_angle = angles.minCoeff(); + auto max_angle = angles.maxCoeff(); + logger().info("Before Min angle: {}, Max angle: {}", min_angle, max_angle); + report["before"]["min_angle"] = min_angle; + report["before"]["max_angle"] = max_angle; + report["before"]["avg_angle"] = angles.mean(); + } + // min/max doublearea + { + Eigen::VectorXd double_areas; + igl::doublearea(inV, inF, double_areas); + auto min_da = double_areas.minCoeff(); + auto max_da = double_areas.maxCoeff(); + logger().info("Before Min double area: {}, Max double area: {}", min_da, max_da); + report["before"]["min_da"] = min_da; + report["before"]["max_da"] = max_da; + report["before"]["avg_da"] = double_areas.mean(); + } + } + + box_minmax = std::pair(inV.colwise().minCoeff(), inV.colwise().maxCoeff()); double diag = (box_minmax.first - box_minmax.second).norm(); - const double envelope_size = env_rel * diag; + if (eps < 0) { + eps = eps_rel * diag; + } igl::Timer timer; + wmtk::logger().info("Total frozen vertices: {}", fixed_vertices.size()); + UniformRemeshing m(verts, num_threads, !sample_envelope); - m.create_mesh(verts.size(), tris, modified_nonmanifold_v, freeze_boundary, envelope_size); + m.set_feature_edges(feature_edge_list); + // m.set_feature_vertices(feature_vertices); + m.create_mesh(verts.size(), tris, fixed_vertices, freeze_boundary, eps); + m.set_patch_ids(face_patch_ids); - { + + if (length_factor < 0) { std::vector properties = m.average_len_valen(); wmtk::logger().info("before remesh properties: {}", properties); if (length_abs < 0 && length_rel < 0) { @@ -129,29 +265,68 @@ int main(int argc, char** argv) } else if (length_abs < 0) { length_abs = diag * length_rel; } + logger().info("absolute target length: {}", length_abs); + m.set_target_edge_length(length_abs); + } else { + logger().info("Use per-patch length with factor {}", length_factor); + m.set_per_patch_target_edge_length(length_factor); + } + auto properties = m.average_len_valen(); + report["before"]["avg_length"] = properties[0]; + report["before"]["max_length"] = properties[1]; + report["before"]["min_length"] = properties[2]; + report["before"]["avg_valence"] = properties[3]; + report["before"]["max_valence"] = properties[4]; + report["before"]["min_valence"] = properties[5]; + + // Write an initial report so downstream code (e.g. write_vtu inside run_remeshing) + // can update/append fields like after min/max angle and double-area. + { + std::ofstream fout(report_file); + fout << std::setw(4) << report; } - logger().info("absolute target length: {}", length_abs); timer.start(); - run_remeshing(input_path, length_abs, output, m, itrs); + run_remeshing(input_path, output, m, itrs, debug_output); timer.stop(); - const std::string report_file = json_params["report"]; - if (!report_file.empty()) { + // Reload report in case downstream code (e.g. write_vtu) has updated it. + try { + std::ifstream fin(report_file); + if (fin) { + fin >> report; + } + } catch (const std::exception& e) { + logger().warn("Failed to reload report file {}: {}", report_file, e.what()); + } + + if (!report.is_object()) { + report = nlohmann::json::object(); + } + if (!report.contains("before") || !report["before"].is_object()) { + report["before"] = nlohmann::json::object(); + } + if (!report.contains("after") || !report["after"].is_object()) { + report["after"] = nlohmann::json::object(); + } + + properties = m.average_len_valen(); + report["after"]["nV"] = static_cast(m.vert_capacity()); + report["after"]["nF"] = static_cast(m.tri_capacity()); + report["after"]["avg_length"] = properties[0]; + report["after"]["max_length"] = properties[1]; + report["after"]["min_length"] = properties[2]; + report["after"]["avg_valence"] = properties[3]; + report["after"]["max_valence"] = properties[4]; + report["after"]["min_valence"] = properties[5]; + + report["time_sec"] = timer.getElapsedTimeInSec(); + + // Persist final merged report. + { std::ofstream fout(report_file); - nlohmann::json report; - std::vector properties = m.average_len_valen(); - report["avg_length"] = properties[0]; - report["max_length"] = properties[1]; - report["min_length"] = properties[2]; - report["avg_valence"] = properties[3]; - report["max_valence"] = properties[4]; - report["min_valence"] = properties[5]; - report["target_length"] = length_abs; - report["time_sec"] = timer.getElapsedTimeInSec(); fout << std::setw(4) << report; - fout.close(); } return 0; -} +} \ No newline at end of file diff --git a/app/remeshing/app/remeshing_spec.hpp b/app/remeshing/app/remeshing_spec.hpp index 6be0fa061a..8df653059b 100644 --- a/app/remeshing/app/remeshing_spec.hpp +++ b/app/remeshing/app/remeshing_spec.hpp @@ -13,12 +13,15 @@ nlohmann::json remeshing_spec = R"( "use_sample_envelope", "num_threads", "max_iterations", + "eps", "eps_rel", "length_rel", "length_abs", + "length_factor", "freeze_boundary", "log_file", - "report" + "report", + "DEBUG_output" ] }, { @@ -50,6 +53,12 @@ nlohmann::json remeshing_spec = R"( "default": 3, "doc": "Maximum iterations before stopping." }, + { + "pointer": "/eps", + "type": "float", + "default": -1, + "doc": "Envelope thickness" + }, { "pointer": "/eps_rel", "type": "float", @@ -66,7 +75,13 @@ nlohmann::json remeshing_spec = R"( "pointer": "/length_abs", "type": "float", "default": -1, - "doc": "Absolute target edge length. If negative, relative length is used to compute the absolute." + "doc": "Absolute target edge length. If negative, relative length is used to compute the absolute. If this is negative as well, the average edge length of the input is used as target." + }, + { + "pointer": "/length_factor", + "type": "float", + "default": -1, + "doc": "Use the per-patch maximum edge length and multiply that with this factor. For coarsening, use a factor >1, for refining <1." }, { "pointer": "/freeze_boundary", @@ -85,6 +100,12 @@ nlohmann::json remeshing_spec = R"( "type": "string", "default": "", "doc": "A JSON file that stores information about the result and the method execution, e.g., runtime." + }, + { + "pointer": "/DEBUG_output", + "type": "bool", + "default": false, + "doc": "Write the mesh as debug_{}.vtu after every operation." } ] )"_json; diff --git a/app/remeshing/app/remeshing_spec.json b/app/remeshing/app/remeshing_spec.json index eca0c6f5a1..e60bcfd288 100644 --- a/app/remeshing/app/remeshing_spec.json +++ b/app/remeshing/app/remeshing_spec.json @@ -8,12 +8,15 @@ "use_sample_envelope", "num_threads", "max_iterations", + "eps", "eps_rel", "length_rel", "length_abs", + "length_factor", "freeze_boundary", "log_file", - "report" + "report", + "DEBUG_output" ] }, { @@ -45,6 +48,12 @@ "default": 3, "doc": "Maximum iterations before stopping." }, + { + "pointer": "/eps", + "type": "float", + "default": -1, + "doc": "Envelope thickness" + }, { "pointer": "/eps_rel", "type": "float", @@ -63,6 +72,12 @@ "default": -1, "doc": "Absolute target edge length. If negative, relative length is used to compute the absolute. If this is negative as well, the average edge length of the input is used as target." }, + { + "pointer": "/length_factor", + "type": "float", + "default": -1, + "doc": "Use the per-patch maximum edge length and multiply that with this factor. For coarsening, use a factor >1, for refining <1." + }, { "pointer": "/freeze_boundary", "type": "bool", @@ -80,5 +95,11 @@ "type": "string", "default": "", "doc": "A JSON file that stores information about the result and the method execution, e.g., runtime." + }, + { + "pointer": "/DEBUG_output", + "type": "bool", + "default": false, + "doc": "Write the mesh as debug_{}.vtu after every operation." } ] diff --git a/app/remeshing/convert_pkl_to_wmtk.py b/app/remeshing/convert_pkl_to_wmtk.py new file mode 100644 index 0000000000..68bbea5127 --- /dev/null +++ b/app/remeshing/convert_pkl_to_wmtk.py @@ -0,0 +1,165 @@ +import igl +import pickle +import numpy as np +import argparse +import json + +# from abs.utils import * # was used to save obj mesh but IGL can do that as well + + +# default paths +fused_path = "fused.pkl" +mesh_path = "mesh.obj" +# feature_edge_vertex_path = "feature_edge_vertex.json" +# feature_vertex_edge_path = "feature_vertex_edge.json" +# corner_path = "corner.json" +patches_path = "features.json" +remesh_json_path = "remesh.json" + +if __name__ == "__main__": + # read arguments + parser = argparse.ArgumentParser() + parser.add_argument("output_dir", type=str, help="output directory") + parser.add_argument( + "--fused", type=str, default=None, help="input fused pkl file" + ) + parser.add_argument( + "--mesh", type=str, default=None, help="output mesh obj file" + ) + # parser.add_argument( + # "--feature_edge_vertex", type=str, default=feature_edge_vertex_path + # ) + # parser.add_argument( + # "--feature_vertex_edge", type=str, default=feature_vertex_edge_path + # ) + # parser.add_argument("--corner", type=str, default=corner_path) + parser.add_argument( + "--patches", type=str, default=None, help="output patches json file" + ) + parser.add_argument( + "--remesh_json", + type=str, + default=None, + help="output remeshing config json file (for -j remesh.json)", + ) + + parser.add_argument("--len_rel", type=float, default=0.05, help="relative length for remeshing") + args = parser.parse_args() + if args.fused: + fused_path = args.fused + print("using fused path: ", args.fused) + else: + fused_path = args.output_dir + "/" + fused_path + + if args.mesh: + mesh_path = args.mesh + else: + mesh_path = args.output_dir + "/" + mesh_path + # feature_edge_vertex_path = args.feature_edge_vertex + # feature_vertex_edge_path = args.feature_vertex_edge + # corner_path = args.corner + if args.patches: + patches_path = args.patches + else: + patches_path = args.output_dir + "/" + patches_path + + if args.remesh_json: + remesh_json_path = args.remesh_json + else: + remesh_json_path = args.output_dir + "/" + remesh_json_path + + + with open(fused_path, "rb") as file: + data = pickle.load(file) + + V = np.array(data[0]) + F = np.array(data[1]) + + # patch: is the face fid: triangle id + patch2fid = data[2] + + # feature edge in the CAD and [list of vertex ids in mesh] + edge2vidschain = data[3] + + # corner is the CAD vertex and the vids is the vertex id in the mesh + corner2vids = data[4] + + # save_obj_mesh(mesh_path, V, F) + igl.write_obj(mesh_path, V, F) + print("writing mesh to ", mesh_path) + + seg2edge = {} + degenerate_feature_edges_to_Cid = {} + for k in edge2vidschain: + if len(edge2vidschain[k]) < 2: + seg2edge[(edge2vidschain[k][0], edge2vidschain[k][0])] = k + degenerate_feature_edges_to_Cid[k] = edge2vidschain[k][0] + continue + for i in range(len(edge2vidschain[k]) - 1): + seg = (edge2vidschain[k][i], edge2vidschain[k][i + 1]) + if seg[1] < seg[0]: + seg = (seg[1], seg[0]) + seg2edge[seg] = k + + vid_to_corner = {} + for c in corner2vids: + vid_to_corner[corner2vids[c]] = c + for k, vid in degenerate_feature_edges_to_Cid.items(): + if vid not in vid_to_corner.keys(): + raise RuntimeError("degenerate feature edge vid not in corner vids") + degenerate_feature_edges_to_Cid[k] = vid_to_corner[vid] + + # with open(feature_edge_vertex_path, "w") as f: + # json.dump(edge2vidschain, f, indent=4) + + seg2edge_tmp = {f"{k[0]},{k[1]}": v for k, v in seg2edge.items()} + + fid2patch = {} + for p in patch2fid: + for fid in patch2fid[p]: + fid2patch[fid] = p + + # with open(feature_vertex_edge_path, "w") as f: + # json.dump(seg2edge_tmp, f, indent=2) + + # with open(corner_path, "w") as f: + # json.dump(corner2vids, f, indent=4) + + corners = {} + for c in corner2vids: + corners[c] = corner2vids[c] + + patches_json = {} + patches_json["fid2patch"] = fid2patch + patches_json["edge2seg"] = seg2edge_tmp + patches_json["degenerate_feature_to_cid"] = degenerate_feature_edges_to_Cid + patches_json["corner_vids"] = corners + + with open(patches_path, "w") as f: + print("writing patche info to ", patches_path) + json.dump(patches_json, f, indent=4) + + # Write a default remeshing config next to the produced mesh/patches so the + # remeshing app can be called with: `-j remesh.json`. + remesh_config = { + "input": mesh_path, + "output": args.output_dir + "/remeshed", + "patches": patches_path, + "log_file": "", + "report": "", + "num_threads": 0, + "max_iterations": 3, + "eps_rel": 0.001, + "use_sample_envelope": True, + "freeze_boundary": False, + "length_factor": -1.0, + "length_abs": -1.0, + "length_rel": args.len_rel, + "DEBUG_output": False, + } + + with open(remesh_json_path, "w") as f: + print("writing remesh config to ", remesh_json_path) + json.dump(remesh_config, f, indent=4) + + print("===== conversion done =====") diff --git a/app/remeshing/remesh_pipeline_one_model.py b/app/remeshing/remesh_pipeline_one_model.py new file mode 100644 index 0000000000..250dfbb578 --- /dev/null +++ b/app/remeshing/remesh_pipeline_one_model.py @@ -0,0 +1,146 @@ +#!/usr/bin/env python3 +import os +import sys +import subprocess +from pathlib import Path + +# ------------------------------------------------- +# Args +# ------------------------------------------------- +if len(sys.argv) != 3 and len(sys.argv) != 4: + print( + "Usage:\n" + " python remesh_one_model.py \n" + ) + sys.exit(1) + +DATASET = sys.argv[1] +MODEL_NAME = sys.argv[2] +if len(sys.argv) >= 4: + PROJ_DIR = sys.argv[3] +else: + if DATASET == "ABC": + PROJ_DIR = "/scratch/yz6675/felix/CAD_2025/" + elif DATASET == "fusion": + PROJ_DIR = "/scratch/yz6675/CAD_2024/" + elif DATASET == "test": + PROJ_DIR= "/Users/yunfanzhou/CAD_2024/" + else: + raise RuntimeError("Please provide project dir") + + +# ------------------------------------------------- +# Paths (adjust once) +# ------------------------------------------------- +if DATASET == "ABC": + RESULT_ROOT = os.path.join( + PROJ_DIR, "result_abc_0", "abc_0000_step_v00" + ) + CAD_INPUT_ROOT = os.path.join( + PROJ_DIR, "abc_0000_step_v00" + ) +elif DATASET == "fusion": + RESULT_ROOT = os.path.join( + PROJ_DIR, "myresult/fusion/a1.0.0_00/" + ) + CAD_INPUT_ROOT = os.path.join( + PROJ_DIR, "data/fusion/a1.0.0_00/" + ) +elif DATASET == "test": + PROJ_ROOT = "/Volumes/myssd/CAD_2024_ssd/" + RESULT_ROOT = os.path.join(PROJ_ROOT, "my_results/abc_test/") + CAD_INPUT_ROOT = os.path.join(PROJ_ROOT, "abc_20_data/") + print("testing mode with hard-coded CAD_INPUT_ROOT and RESULT_ROOT") +else: + raise RuntimeError("Unknown dataset: " + DATASET) +MODEL_RESULT_DIR = os.path.join(RESULT_ROOT, MODEL_NAME) + +FUSED_PKL = os.path.join(MODEL_RESULT_DIR, "fused.pkl") +REMESH_JSON = os.path.join(MODEL_RESULT_DIR, "remesh.json") +REMESH_DONE = os.path.join(MODEL_RESULT_DIR, "remesh.done") +CONVERT_SCRIPT = os.path.join( + PROJ_DIR, + "remeshing", + "wildmeshing-toolkit", + "app", + "remeshing", + "convert_pkl_to_wmtk.py", +) + +REMESH_APP = os.path.join( + PROJ_DIR, + "remeshing", + "wildmeshing-toolkit", + "build_stable", + "app", + "remeshing_app", +) + +POST_PROCESS_APP = os.path.join( + PROJ_DIR, + "cad_meshing", + "metaprogramming", + "convert_remeshed_vtus.py", +) + +# ------------------------------------------------- +def run(cmd, cwd=None): + print(">>>", " ".join(cmd)) + subprocess.check_call(cmd, cwd=cwd) + +# ------------------------------------------------- +print(f"=== Remeshing {MODEL_NAME} ===") + +# ------------------------- +# Stage 0: Sanity checks +# ------------------------- +if not os.path.isdir(MODEL_RESULT_DIR): + raise RuntimeError(f"Missing model result dir: {MODEL_RESULT_DIR}") + +if not os.path.exists(FUSED_PKL): + print("[Skip] fused.pkl not found") + sys.exit(0) + +# ------------------------- +# Stage 1: fused.pkl -> remesh.json +# ------------------------- +if not os.path.exists(REMESH_JSON): + run([ + "python", + CONVERT_SCRIPT, + MODEL_RESULT_DIR, + ]) +else: + print("[Stage 1] remesh.json exists, skipped") + +if not os.path.exists(REMESH_JSON): + raise RuntimeError("Failed to generate remesh.json") + +# ------------------------- +# Stage 2: remeshing_app +# ------------------------- +if not os.path.exists(REMESH_DONE): + run([ + REMESH_APP, + "-j", + REMESH_JSON, + ]) + + # success marker + Path(REMESH_DONE).touch() +else: + print("[Stage 2] remeshing already done, skipped") + +print(f"=== Done remeshing {MODEL_NAME} ===") + +# ------------------------- +# Stage 3: verify and render +# ------------------------- +if os.path.exists(REMESH_DONE): + run(["python", + POST_PROCESS_APP, + MODEL_NAME, + CAD_INPUT_ROOT, + RESULT_ROOT]) + +print(f"=== Done remeshing post-processing {MODEL_NAME} ===") \ No newline at end of file diff --git a/app/remeshing/src/remeshing/UniformRemeshing.cpp b/app/remeshing/src/remeshing/UniformRemeshing.cpp index 525e5d21e3..b44d6a87f3 100644 --- a/app/remeshing/src/remeshing/UniformRemeshing.cpp +++ b/app/remeshing/src/remeshing/UniformRemeshing.cpp @@ -7,11 +7,14 @@ #include #include #include +#include #include +#include #include -using namespace app::remeshing; + using namespace wmtk; +namespace { auto renew = [](auto& m, auto op, auto& tris) { auto edges = m.new_edges_after(tris); auto optup = std::vector>(); @@ -25,6 +28,12 @@ auto edge_locker = [](auto& m, const auto& e, int task_id) { return m.try_set_edge_mutex_two_ring(e, task_id); }; +static int debug_print_counter = 0; +} // namespace + +namespace app::remeshing { + + UniformRemeshing::UniformRemeshing( std::vector _m_vertex_positions, int num_threads, @@ -34,6 +43,8 @@ UniformRemeshing::UniformRemeshing( m_envelope.use_exact = use_exact; p_vertex_attrs = &vertex_attrs; + p_edge_attrs = &edge_attrs; + p_face_attrs = &face_attrs; vertex_attrs.resize(_m_vertex_positions.size()); @@ -44,53 +55,153 @@ UniformRemeshing::UniformRemeshing( void UniformRemeshing::create_mesh( size_t n_vertices, const std::vector>& tris, - const std::vector& frozen_verts, + const std::vector>& frozen_verts, bool m_freeze, double eps) { wmtk::TriMesh::init(n_vertices, tris); - std::vector V(n_vertices); - std::vector F(tris.size()); - for (auto i = 0; i < V.size(); i++) { - V[i] = vertex_attrs[i].pos; - } - for (int i = 0; i < F.size(); ++i) F[i] << tris[i][0], tris[i][1], tris[i][2]; - if (eps > 0) { - m_envelope.init(V, F, eps); - m_has_envelope = true; - } else - m_envelope.init(V, F, 0.0); // TODO: this should not be here partition_mesh_morton(); + for (const auto& [v, corner_id] : frozen_verts) { + vertex_attrs[v].corner_id = corner_id; + } + if (m_freeze) { - for (auto v : frozen_verts) { - vertex_attrs[v].freeze = true; - } - for (auto e : get_edges()) { + for (const Tuple& e : get_edges()) { if (is_boundary_edge(e)) { - vertex_attrs[e.vid(*this)].freeze = true; - vertex_attrs[e.switch_vertex(*this).vid(*this)].freeze = true; + vertex_attrs[e.vid(*this)].corner_id = 100000; + vertex_attrs[e.switch_vertex(*this).vid(*this)].corner_id = 100000; } } } + edge_attrs.resize(tri_capacity() * 3); + initialize_feature_edges(); + + // init envelopes + { + std::vector V(n_vertices); + std::vector F(tris.size()); + const auto feature_edges = + get_edges_by_condition([](const EdgeAttributes& e) { return e.is_feature; }); + std::vector E(feature_edges.size()); + + for (int i = 0; i < V.size(); i++) { + V[i] = vertex_attrs[i].pos; + } + for (int i = 0; i < F.size(); ++i) { + F[i] = Vector3i(tris[i][0], tris[i][1], tris[i][2]); + } + for (int i = 0; i < E.size(); ++i) { + E[i] = Vector2i(feature_edges[i][0], feature_edges[i][1]); + } + if (eps > 0) { + m_envelope.init(V, F, eps); + if (!feature_edges.empty()) { + m_feature_envelope.init(V, E, eps); + } + m_has_envelope = true; + } else { + m_envelope.init(V, F, 0.0); + if (!feature_edges.empty()) { + m_feature_envelope.init(V, E, 0.0); + } + } + } +} + +bool UniformRemeshing::cache_edge_positions(const Tuple& t) +{ + auto& cache = collapse_info_cache.local(); + cache.link_vertices.clear(); + + const size_t v0 = t.vid(*this); + const size_t v1 = t.switch_vertex(*this).vid(*this); + + cache.v0p = vertex_attrs.at(v0).pos; + cache.v1p = vertex_attrs.at(v1).pos; + cache.v0_tal = vertex_attrs.at(v0).tal; + cache.v1_tal = vertex_attrs.at(v1).tal; + cache.v0_corner_id = vertex_attrs.at(v0).corner_id; + cache.v1_corner_id = vertex_attrs.at(v1).corner_id; + cache.partition_id = vertex_attrs.at(v0).partition_id; + + cache.v0 = v0; + cache.v1 = v1; + cache.is_feature_edge = edge_attrs.at(t.eid(*this)).is_feature; + + auto add_link_vertices = [this, &v0, &v1, &cache](const Tuple& t_) { + // use a dummy to make sure the default value for is_feature is used + EdgeAttributes dummy; + + const size_t v2 = t_.switch_edge(*this).switch_vertex(*this).vid(*this); + const size_t e0 = eid_from_vids(v0, v2); + const size_t e1 = eid_from_vids(v1, v2); + const int e0_is_feature = edge_attrs.at(e0).is_feature; + const int e1_is_feature = edge_attrs.at(e1).is_feature; + if (e0_is_feature != dummy.is_feature && e1_is_feature != dummy.is_feature) { + return false; // cannot merge two feature edges + } + if (e0_is_feature != dummy.is_feature) { + cache.link_vertices.emplace_back(v2, e0_is_feature); + } else if (e1_is_feature != dummy.is_feature) { + cache.link_vertices.emplace_back(v2, e1_is_feature); + } else { + cache.link_vertices.emplace_back(v2, dummy.is_feature); + } + return true; + }; + + if (!add_link_vertices(t)) { + return false; + } + for (const Tuple& t_opp : t.switch_faces(*this)) { + if (!add_link_vertices(t_opp)) { + return false; + } + } + + return true; } -void UniformRemeshing::cache_edge_positions(const Tuple& t) +std::vector> UniformRemeshing::get_edges_by_condition( + std::function cond) const { - position_cache.local().v1p = vertex_attrs[t.vid(*this)].pos; - position_cache.local().v2p = vertex_attrs[t.switch_vertex(*this).vid(*this)].pos; - position_cache.local().partition_id = vertex_attrs[t.vid(*this)].partition_id; + std::vector> res; + for (const Tuple& e : get_edges()) { + auto eid = e.eid(*this); + if (cond(edge_attrs[eid])) { + size_t v0 = e.vid(*this); + size_t v1 = e.switch_vertex(*this).vid(*this); + res.emplace_back(std::array{{v0, v1}}); + } + } + return res; } bool UniformRemeshing::invariants(const std::vector& new_tris) { + for (const Tuple& t : new_tris) { + const double q = get_quality(t); + if (q == 0) { + return false; + } + } + if (m_has_envelope) { for (auto& t : new_tris) { std::array tris; auto vs = oriented_tri_vertices(t); - for (auto j = 0; j < 3; j++) tris[j] = vertex_attrs[vs[j].vid(*this)].pos; + + for (auto j = 0; j < 3; j++) { + const auto vid = vs[j].vid(*this); + tris[j] = vertex_attrs[vid].pos; + + if (vertex_attrs[vid].is_feature && m_feature_envelope.is_outside(tris[j])) { + return false; + } + } if (m_envelope.is_outside(tris)) { return false; } @@ -215,21 +326,120 @@ std::vector UniformRemeshing::new_edges_after( bool UniformRemeshing::swap_edge_before(const Tuple& t) { + // Do not swap feature edges + if (is_feature_edge(t)) { + return false; + } if (!TriMesh::swap_edge_before(t)) return false; - if (vertex_attrs[t.vid(*this)].freeze && vertex_attrs[t.switch_vertex(*this).vid(*this)].freeze) + + const size_t v0 = t.vid(*this); + const size_t v1 = t.switch_vertex(*this).vid(*this); + + // do not swap if vertex has valence 3 + if (get_valence_for_vertex(v0) == 3 || get_valence_for_vertex(v1) == 3) { + return false; + } + + // if (vertex_attrs[v0].is_freeze && vertex_attrs[v1].is_freeze) { + // return false; + // } + + auto& cache = swap_info_cache.local(); + cache.ring_edge_attrs.clear(); + + cache.v0 = v0; + cache.v1 = v1; + + cache.v0p = vertex_attrs.at(cache.v0).pos; + cache.v1p = vertex_attrs.at(cache.v1).pos; + + cache.is_feature_edge = 0; // must be false, was checked before + + cache.v2 = t.switch_edge(*this).switch_vertex(*this).vid(*this); + + size_t face_count = 0; + for (const Tuple& t_opp : t.switch_faces(*this)) { + cache.v3 = t_opp.switch_edge(*this).switch_vertex(*this).vid(*this); + ++face_count; + if (face_count > 1) { + // Non-manifold edge (????) + return false; + } + } + if (face_count == 0) { + // Boundary edge (???) return false; + } + + // normal check + { + const size_t v2 = cache.v2; + const size_t v3 = cache.v3; + const Vector3d& p0 = cache.v0p; + const Vector3d& p1 = cache.v1p; + const Vector3d& p2 = vertex_attrs.at(v2).pos; + const Vector3d& p3 = vertex_attrs.at(v3).pos; + + // current + Vector3d n0 = (p1 - p0).cross(p2 - p0); + Vector3d n1 = (p3 - p0).cross(p1 - p0); + if (n0.norm() > 1e-12 && n1.norm() > 1e-12) { + n0.normalize(); + n1.normalize(); + const double a_before = n0.dot(n1); + // after swap + Vector3d n2 = (p2 - p3).cross(p0 - p3).normalized(); + Vector3d n3 = (p1 - p3).cross(p2 - p3).normalized(); + const double a_after = n2.dot(n3); + + if (a_after < a_before && a_after < -0.5) { + return false; + } + } + } + + const std::array ring_edges = { + simplex::Edge(cache.v0, cache.v2), + simplex::Edge(cache.v1, cache.v2), + simplex::Edge(cache.v0, cache.v3), + simplex::Edge(cache.v1, cache.v3)}; + + for (const auto& e : ring_edges) { + const auto& vs = e.vertices(); + const size_t eid = eid_from_vids(vs[0], vs[1]); + cache.ring_edge_attrs[e] = edge_attrs.at(eid); + } return true; } +// bool UniformRemeshing::swap_edge_after(const TriMesh::Tuple& t) +// { +// const auto& cache = swap_info_cache.local(); + +// std::vector tris; +// tris.push_back(t); +// tris.push_back(t.switch_edge(*this)); +// return true; +// } + bool UniformRemeshing::swap_edge_after(const TriMesh::Tuple& t) { - std::vector tris; - tris.push_back(t); - tris.push_back(t.switch_edge(*this)); + const auto& cache = swap_info_cache.local(); + + for (const auto& [edge, attrs] : cache.ring_edge_attrs) { + const auto& vs = edge.vertices(); + const size_t eid = eid_from_vids(vs[0], vs[1]); + edge_attrs[eid] = attrs; + } + const size_t e_new = eid_from_vids(cache.v2, cache.v3); + + edge_attrs[e_new].is_feature = 0; + return true; } + std::vector UniformRemeshing::replace_edges_after_split( const std::vector& tris, const size_t vid_threshold) const @@ -238,7 +448,7 @@ std::vector UniformRemeshing::replace_edges_after_split( // only push back those edges which are already present, but invalidated by hash mechanism. std::vector new_edges; new_edges.reserve(tris.size()); - for (auto t : tris) { + for (const Tuple& t : tris) { auto tmptup = (t.switch_vertex(*this)).switch_edge(*this); if (tmptup.vid(*this) < vid_threshold && (tmptup.switch_vertex(*this)).vid(*this) < vid_threshold) @@ -266,73 +476,308 @@ std::vector UniformRemeshing::new_sub_edges_after_split( bool UniformRemeshing::collapse_edge_before(const Tuple& t) { - if (!TriMesh::collapse_edge_before(t)) return false; - if (vertex_attrs[t.vid(*this)].freeze || vertex_attrs[t.switch_vertex(*this).vid(*this)].freeze) + size_t v0 = t.vid(*this); + size_t v1 = t.switch_vertex(*this).vid(*this); + + // if (is_feature_edge(t)) { + // return false; + // } + + if (!is_feature_edge(t) && (vertex_attrs[v0].is_feature || vertex_attrs[v1].is_feature)) { return false; - cache_edge_positions(t); + } + + // if (is_feature_edge(t) && !vertex_attrs[v0].is_feature && vertex_attrs[v1].is_feature) + // return false; + + + // if (vertex_attrs[v0].corner_id >= 0 || vertex_attrs[v1].corner_id >= 0) { + // return false; + // } + if (vertex_attrs[v0].corner_id >= 0 && vertex_attrs[v1].corner_id >= 0) { + return false; + } + + if (!TriMesh::collapse_edge_before(t)) { + return false; + } + if (!cache_edge_positions(t)) { + return false; + } + + //// for debug + // if (v0 == 35465 && v1 == 36590) { + // logger().error("DEBUG ({},{})", v0, v1); + // logger().error("v0.is_feature = {}", vertex_attrs[v0].is_feature); + // logger().error("v1.is_feature = {}", vertex_attrs[v1].is_feature); + // logger().error("v0.corner_id = {}", vertex_attrs[v0].corner_id); + // logger().error("v1.corner_id = {}", vertex_attrs[v1].corner_id); + // logger().error("is_feature_edge = {}", is_feature_edge(t)); + // logger().error("t.is_feature = {}", edge_attrs[t.eid(*this)].is_feature); + // + // { // DEBUG + // logger().error("DEBUG v0 = {}", v0); + // const Tuple v = tuple_from_vertex(v0); + // for (const Tuple& e : get_one_ring_edges_for_vertex(v)) { + // simplex::Edge s = simplex_from_edge(e); + // logger().error("({}) = {}", s.vertices(), edge_attrs[e.eid(*this)].is_feature); + // } + // } + // { // DEBUG + // logger().error("DEBUG v1 = {}", v1); + // const Tuple v = tuple_from_vertex(v1); + // for (const Tuple& e : get_one_ring_edges_for_vertex(v)) { + // simplex::Edge s = simplex_from_edge(e); + // logger().error("({}) = {}", s.vertices(), edge_attrs[e.eid(*this)].is_feature); + // } + // } + //} + return true; } bool UniformRemeshing::collapse_edge_after(const TriMesh::Tuple& t) { - const Eigen::Vector3d p = (position_cache.local().v1p + position_cache.local().v2p) / 2.0; - auto vid = t.vid(*this); - vertex_attrs[vid].pos = p; - vertex_attrs[vid].partition_id = position_cache.local().partition_id; + const auto& cache = collapse_info_cache.local(); + + auto& attr = vertex_attrs[t.vid(*this)]; + + if (cache.v0_corner_id >= 0) { + attr.pos = cache.v0p; + attr.corner_id = cache.v0_corner_id; + } else if (cache.v1_corner_id >= 0) { + attr.pos = cache.v1p; + attr.corner_id = cache.v1_corner_id; + } else if (cache.is_feature_edge || !vertex_attrs.at(cache.v1).is_feature) { + // collase to midpoint if v1 is not a feature + const Eigen::Vector3d p = (cache.v0p + cache.v1p) / 2.0; + attr.pos = p; + } + attr.partition_id = cache.partition_id; + + + /* + * I am not sure if it is better to keep the target edge length of the remaining vertex or to + * use the min of the two. The min seems to be quite aggressive. + */ + // attr.tal = std::min(cache.v0_tal, cache.v1_tal); + + // update the edge attributes of the surrounding edges + for (const auto [v2, f] : cache.link_vertices) { + const size_t eid = eid_from_vids(t.vid(*this), v2); + edge_attrs[eid].is_feature = f; + } + + // sanity check (for debugging) + if (vertex_attrs.at(cache.v1).is_feature && vertex_attrs.at(cache.v1).corner_id < 0) { + // const Tuple v = tuple_from_vertex(55432); + size_t counter = 0; + for (const Tuple& e : get_one_ring_edges_for_vertex(t)) { + if (edge_attrs[e.eid(*this)].is_feature) { + counter++; + } + } + if (counter != 2) { + logger().error( + "v1 = {}, corner_id = {}", + cache.v1, + vertex_attrs.at(cache.v1).corner_id); + logger().error("Vertex now endpoint, after collapsing ({},{})", cache.v0, cache.v1); + for (const Tuple& e : get_one_ring_edges_for_vertex(t)) { + simplex::Edge s = simplex_from_edge(e); + logger().error("({}) = {}", s.vertices(), edge_attrs[e.eid(*this)].is_feature); + } + log_and_throw_error("Sanity check for feature edges failed!"); + } + } return true; } + bool UniformRemeshing::split_edge_before(const Tuple& t) { - if (!TriMesh::split_edge_before(t)) return false; - cache_edge_positions(t); + if (!TriMesh::split_edge_before(t)) { + return false; + } + + auto& cache = split_info_cache.local(); + cache.edge_attrs.clear(); + cache.face_attrs.clear(); + + cache.v0 = t.vid(*this); + cache.v1 = t.switch_vertex(*this).vid(*this); + + cache.v0p = vertex_attrs.at(cache.v0).pos; + cache.v1p = vertex_attrs.at(cache.v1).pos; + cache.v0_tal = vertex_attrs.at(cache.v0).tal; + cache.v1_tal = vertex_attrs.at(cache.v1).tal; + cache.partition_id = vertex_attrs.at(cache.v0).partition_id; + cache.is_feature_edge = edge_attrs[t.eid(*this)].is_feature; + + // gather all surrounding edges + std::vector edges; + edges.reserve(4); + { + const size_t v_opp = t.switch_edge(*this).switch_vertex(*this).vid(*this); + edges.emplace_back(cache.v0, v_opp); + edges.emplace_back(cache.v1, v_opp); + cache.face_attrs[v_opp] = face_attrs.at(t.fid(*this)); + } + for (const Tuple& t_opp : t.switch_faces(*this)) { + const size_t v_opp = t_opp.switch_edge(*this).switch_vertex(*this).vid(*this); + edges.emplace_back(cache.v0, v_opp); + edges.emplace_back(cache.v1, v_opp); + cache.face_attrs[v_opp] = face_attrs.at(t_opp.fid(*this)); + } + + // save all edge attributes + for (const auto& e : edges) { + const size_t eid = eid_from_vids(e.vertices()[0], e.vertices()[1]); + cache.edge_attrs[e] = edge_attrs.at(eid); + } + return true; } - bool UniformRemeshing::split_edge_after(const TriMesh::Tuple& t) { - const Eigen::Vector3d p = (position_cache.local().v1p + position_cache.local().v2p) / 2.0; - auto vid = t.switch_vertex(*this).vid(*this); - vertex_attrs[vid].pos = p; - vertex_attrs[vid].partition_id = position_cache.local().partition_id; + const auto& cache = split_info_cache.local(); + + const Eigen::Vector3d p = 0.5 * (cache.v0p + cache.v1p); + // Eigen::Vector3d after_project; + // if (cache.is_feature_edge) { + // m_feature_envelope.nearest_point(p, after_project); + // } else { + // m_envelope.nearest_point(p, after_project); + // } + + const size_t vnew = t.switch_vertex(*this).vid(*this); + auto& vnew_attrs = vertex_attrs[vnew]; + vnew_attrs.pos = p; + vnew_attrs.partition_id = cache.partition_id; + vnew_attrs.is_feature = cache.is_feature_edge; + vnew_attrs.tal = 0.5 * (cache.v0_tal + cache.v1_tal); + + // update the edge attributes of the surrounding edges + std::vector vids; + vids.reserve(cache.edge_attrs.size() * 2); + for (const auto& [edge, attrs] : cache.edge_attrs) { + const size_t eid = eid_from_vids(edge.vertices()[0], edge.vertices()[1]); + edge_attrs[eid] = attrs; + vids.push_back(edge.vertices()[0]); + vids.push_back(edge.vertices()[1]); + } + wmtk::vector_unique(vids); + + + // set attributes of new edges + for (size_t vid : vids) { + const size_t e = eid_from_vids(vnew, vid); + if (cache.is_feature_edge && (vid == cache.v0 || vid == cache.v1)) { + edge_attrs[e].is_feature = cache.is_feature_edge; + } else { + edge_attrs[e].is_feature = 0; + } + } + + // set attributes of new faces + for (const auto& [v_opp, attrs] : cache.face_attrs) { + const Tuple f0 = tuple_from_vids(v_opp, vnew, cache.v0); + const Tuple f1 = tuple_from_vids(v_opp, vnew, cache.v1); + face_attrs[f0.fid(*this)] = attrs; + face_attrs[f1.fid(*this)] = attrs; + } + return true; } + bool UniformRemeshing::smooth_before(const Tuple& t) { - if (vertex_attrs[t.vid(*this)].freeze) return false; + if (vertex_attrs[t.vid(*this)].corner_id >= 0) return false; + if (vertex_attrs[t.vid(*this)].is_feature) return false; + return true; } bool UniformRemeshing::smooth_after(const TriMesh::Tuple& t) { auto one_ring_tris = get_one_ring_tris_for_vertex(t); - if (one_ring_tris.size() < 2) { + if (one_ring_tris.size() <= 2) { return false; } + + // store min quality before + double min_q = std::numeric_limits::max(); + for (const Tuple& n : one_ring_tris) { + double q = get_quality(n); + min_q = std::min(min_q, q); + } + Eigen::Vector3d after_smooth = tangential_smooth(t); if (after_smooth.hasNaN()) return false; - vertex_attrs[t.vid(*this)].pos = after_smooth; + // todo: add envelope projection and check here + // if (m_has_envelope) { + // std::array tri; + // for (auto& tri_tup : one_ring_tris) { + // auto vs = oriented_tri_vertices(tri_tup); + // for (auto j = 0; j < 3; j++) + // tri[j] = vertex_attrs[vs[j].vid(*this)].pos; + // if (m_envelope.is_outside(tri)) { + // after_smooth = m_envelope.project_to_envelope(after_smooth); + // break; + // } + // } + + Eigen::Vector3d after_project; + if (vertex_attrs[t.vid(*this)].is_feature) { + m_feature_envelope.nearest_point(after_smooth, after_project); + } else { + m_envelope.nearest_point(after_smooth, after_project); + } + // logger().warn("{} --> {}", after_smooth.transpose(), after_project.transpose()); + + vertex_attrs[t.vid(*this)].pos = after_project; + + // reject if quality after is worse + for (const Tuple& n : one_ring_tris) { + double q = get_quality(n); + if (q < min_q && q < 0.1) { + return false; + } + } + return true; } -double UniformRemeshing::compute_edge_cost_collapse(const TriMesh::Tuple& t, double L) const +double UniformRemeshing::compute_edge_cost_collapse(const TriMesh::Tuple& t) const { - double l = - (vertex_attrs[t.vid(*this)].pos - vertex_attrs[t.switch_vertex(*this).vid(*this)].pos) - .norm(); - if (l < (4. / 5.) * L) return ((4. / 5.) * L - l); + const size_t v0 = t.vid(*this); + const size_t v1 = t.switch_vertex(*this).vid(*this); + + const double tal = 0.5 * (vertex_attrs.at(v0).tal + vertex_attrs.at(v0).tal); + + double l = (vertex_attrs.at(v0).pos - vertex_attrs.at(v1).pos).norm(); + + if (l < (4. / 5.) * tal) { + return ((4. / 5.) * tal - l); + } return -1; } -double UniformRemeshing::compute_edge_cost_split(const TriMesh::Tuple& t, double L) const +double UniformRemeshing::compute_edge_cost_split(const TriMesh::Tuple& t) const { - double l = - (vertex_attrs[t.vid(*this)].pos - vertex_attrs[t.switch_vertex(*this).vid(*this)].pos) - .norm(); - if (l > (4. / 3.) * L) return (l - (4. / 3.) * L); + const size_t v0 = t.vid(*this); + const size_t v1 = t.switch_vertex(*this).vid(*this); + + const double tal = 0.5 * (vertex_attrs.at(v0).tal + vertex_attrs.at(v0).tal); + + double l = (vertex_attrs.at(v0).pos - vertex_attrs.at(v1).pos).norm(); + + if (l > (4. / 3.) * tal) { + return (l - (4. / 3.) * tal); + } return -1; } @@ -345,12 +790,14 @@ double UniformRemeshing::compute_vertex_valence(const TriMesh::Tuple& t) const auto t3 = (t.switch_edge(*this)).switch_vertex(*this); valences[2] = std::make_pair(t3, get_valence_for_vertex(t3)); - if ((t.switch_face(*this)).has_value()) { - auto t4 = (((t.switch_face(*this)).value()).switch_edge(*this)).switch_vertex(*this); + const auto t_opps = t.switch_faces(*this); + + if (t_opps.size() == 1) { + auto t4 = t_opps[0].switch_edge(*this).switch_vertex(*this); valences.emplace_back(t4, get_valence_for_vertex(t4)); } - double cost_before_swap = 0.0; - double cost_after_swap = 0.0; + int cost_before_swap = 0; + int cost_after_swap = 0; // check if it's internal vertex or bondary vertex // navigating starting one edge and getting back to the start @@ -358,21 +805,44 @@ double UniformRemeshing::compute_vertex_valence(const TriMesh::Tuple& t) const for (int i = 0; i < valences.size(); i++) { TriMesh::Tuple vert = valences[i].first; int val = 6; - auto one_ring_edges = get_one_ring_edges_for_vertex(vert); - for (auto edge : one_ring_edges) { + const auto one_ring_edges = get_one_ring_edges_for_vertex(vert); + for (const Tuple& edge : one_ring_edges) { if (is_boundary_edge(edge)) { val = 4; break; } } - cost_before_swap += (double)(valences[i].second - val) * (valences[i].second - val); - cost_after_swap += - (i < 2) ? (double)(valences[i].second - 1 - val) * (valences[i].second - 1 - val) - : (double)(valences[i].second + 1 - val) * (valences[i].second + 1 - val); + cost_before_swap += (valences[i].second - val) * (valences[i].second - val); + cost_after_swap += (i < 2) + ? (valences[i].second - 1 - val) * (valences[i].second - 1 - val) + : (valences[i].second + 1 - val) * (valences[i].second + 1 - val); } - return (cost_before_swap - cost_after_swap); + return (double)(cost_before_swap - cost_after_swap); } +double UniformRemeshing::compute_swap_energy(const Tuple& t) const +{ + const auto t_opps = t.switch_faces(*this); + if (t_opps.size() != 1) { + return -1e50; // don't swap this + } + const size_t v0 = t.vid(*this); + const size_t v1 = t.switch_vertex(*this).vid(*this); + const size_t v2 = t.switch_edge(*this).switch_vertex(*this).vid(*this); + const size_t v3 = t_opps[0].switch_edge(*this).switch_vertex(*this).vid(*this); + + // old: (0,1,2) (0,1,3) + // new: (0,2,3) (1,2,3) + + const double q012 = get_quality({{v0, v1, v2}}); // old + const double q013 = get_quality({{v0, v1, v3}}); // old + const double q023 = get_quality({{v0, v2, v3}}); // new + const double q123 = get_quality({{v1, v2, v3}}); // new + const double before_swap = 1.0 / std::min(q012, q013); + const double after_swap = 1.0 / std::min(q023, q123); + + return before_swap - after_swap; +} std::vector UniformRemeshing::average_len_valen() { double average_len = 0.0; @@ -415,7 +885,7 @@ std::vector UniformRemeshing::new_edges_after_swap(const TriMesh return new_edges; } -bool UniformRemeshing::collapse_remeshing(double L) +bool UniformRemeshing::collapse_remeshing() { igl::Timer timer; timer.start(); @@ -424,7 +894,9 @@ bool UniformRemeshing::collapse_remeshing(double L) for_each_edge([&](auto& tup) { collect_tuples.emplace_back(tup); }); collect_all_ops.reserve(collect_tuples.size()); - for (auto& t : collect_tuples) collect_all_ops.emplace_back("edge_collapse", t); + for (auto& t : collect_tuples) { + collect_all_ops.emplace_back("edge_collapse", t); + } wmtk::logger().info( "***** collapse get edges time *****: {} ms", timer.getElapsedTimeInMilliSec()); @@ -433,13 +905,15 @@ bool UniformRemeshing::collapse_remeshing(double L) auto setup_and_execute = [&](auto& executor) { executor.renew_neighbor_tuples = renew; executor.priority = [&](auto& m, auto _, auto& e) { - return m.compute_edge_cost_collapse(e, L); + return m.compute_edge_cost_collapse(e); }; executor.num_threads = NUM_THREADS; executor.should_renew = [](auto val) { return (val > 0); }; executor.is_weight_up_to_date = [](auto& m, auto& ele) { auto& [val, op, e] = ele; - if (val < 0) return false; // priority is negated. + if (val < 0) { + return false; // priority is negated. + } return true; }; executor(*this, collect_all_ops); @@ -455,7 +929,7 @@ bool UniformRemeshing::collapse_remeshing(double L) return true; } -bool UniformRemeshing::split_remeshing(double L) +bool UniformRemeshing::split_remeshing() { igl::Timer timer; timer.start(); @@ -466,41 +940,51 @@ bool UniformRemeshing::split_remeshing(double L) for_each_edge([&](auto& tup) { collect_tuples.emplace_back(tup); }); collect_all_ops.reserve(collect_tuples.size()); - for (auto& t : collect_tuples) collect_all_ops.emplace_back("edge_split", t); + for (auto& t : collect_tuples) { + collect_all_ops.emplace_back("edge_split", t); + } wmtk::logger().info( "***** split get edges time *****: {} ms", timer.getElapsedTimeInMilliSec()); wmtk::logger().info("size for edges to be split is {}", collect_all_ops.size()); - auto edges2 = tbb::concurrent_vector>(); + // auto edges2 = tbb::concurrent_vector>(); auto setup_and_execute = [&](auto& executor) { vid_threshold = vert_capacity(); executor.num_threads = NUM_THREADS; executor.renew_neighbor_tuples = [&](auto& m, auto op, auto& tris) { count_success++; auto edges = m.replace_edges_after_split(tris, vid_threshold); - for (auto e2 : m.new_sub_edges_after_split(tris)) edges2.emplace_back(op, e2); + // for (const Tuple& e2 : m.new_sub_edges_after_split(tris)) { + // edges2.emplace_back(op, e2); + // } auto optup = std::vector>(); - for (auto& e : edges) optup.emplace_back(op, e); + for (const Tuple& e : edges) { + optup.emplace_back(op, e); + } return optup; }; - executor.priority = [&](auto& m, auto _, auto& e) { - return m.compute_edge_cost_split(e, L); - }; + executor.priority = [&](auto& m, auto _, auto& e) { return m.compute_edge_cost_split(e); }; executor.should_renew = [](auto val) { return (val > 0); }; executor.is_weight_up_to_date = [](auto& m, auto& ele) { auto& [val, op, e] = ele; - if (val < 0) return false; + if (val < 0) { + return false; + } return true; }; // Execute!! - do { - count_success.store(0, std::memory_order_release); - executor(*this, collect_all_ops); - collect_all_ops.clear(); - for (auto& item : edges2) collect_all_ops.emplace_back(item); - edges2.clear(); - } while (count_success.load(std::memory_order_acquire) > 0); + count_success.store(0, std::memory_order_release); + executor(*this, collect_all_ops); + // do { + // count_success.store(0, std::memory_order_release); + // executor(*this, collect_all_ops); + // collect_all_ops.clear(); + // for (auto& item : edges2) { + // collect_all_ops.emplace_back(item); + // } + // edges2.clear(); + // } while (count_success.load(std::memory_order_acquire) > 0); }; if (NUM_THREADS > 0) { @@ -518,7 +1002,9 @@ bool UniformRemeshing::split_remeshing(double L) bool UniformRemeshing::smooth_all_vertices() { auto collect_all_ops = std::vector>(); - for (auto& loc : get_edges()) collect_all_ops.emplace_back("vertex_smooth", loc); + for (auto& loc : get_vertices()) { + collect_all_ops.emplace_back("vertex_smooth", loc); + } auto setup_and_execute = [&](auto& executor) { executor.num_threads = NUM_THREADS; @@ -557,12 +1043,12 @@ bool UniformRemeshing::swap_remeshing() executor.renew_neighbor_tuples = renew; executor.num_threads = NUM_THREADS; executor.priority = [](auto& m, auto op, const Tuple& e) { - return m.compute_vertex_valence(e); + return m.compute_swap_energy(e); }; executor.should_renew = [](auto val) { return (val > 0); }; executor.is_weight_up_to_date = [](auto& m, auto& ele) { auto& [val, _, e] = ele; - auto val_energy = (m.compute_vertex_valence(e)); + auto val_energy = (m.compute_swap_energy(e)); return (val_energy > 1e-5) && ((val_energy - val) * (val_energy - val) < 1e-8); }; executor(*this, collect_all_ops); @@ -622,65 +1108,85 @@ Eigen::Vector3d UniformRemeshing::tangential_smooth(const Tuple& t) if (one_ring_tris.size() < 2) return vertex_attrs[t.vid(*this)].pos; Eigen::Vector3d after_smooth = smooth(t); // get normal and area of each face - auto area = [](auto& m, auto& verts) { - return ((m.vertex_attrs[verts[0].vid(m)].pos - m.vertex_attrs[verts[2].vid(m)].pos) - .cross( - m.vertex_attrs[verts[1].vid(m)].pos - m.vertex_attrs[verts[2].vid(m)].pos)) - .norm() / - 2.0; - }; - auto normal = [](auto& m, auto& verts) { - return ((m.vertex_attrs[verts[0].vid(m)].pos - m.vertex_attrs[verts[2].vid(m)].pos) - .cross( - m.vertex_attrs[verts[1].vid(m)].pos - m.vertex_attrs[verts[2].vid(m)].pos)) - .normalized(); - }; - auto w0 = 0.0; - Eigen::Vector3d n0(0.0, 0.0, 0.0); - for (auto& e : one_ring_tris) { - auto verts = oriented_tri_vertices(e); - w0 += area(*this, verts); - n0 += area(*this, verts) * normal(*this, verts); - } - n0 /= w0; - if (n0.norm() < 1e-10) return vertex_attrs[t.vid(*this)].pos; - n0 = n0.normalized(); - after_smooth += n0 * n0.transpose() * (vertex_attrs[t.vid(*this)].pos - after_smooth); - assert(check_mesh_connectivity_validity()); + // auto area = [](auto& m, auto& verts) { + // return ((m.vertex_attrs[verts[0].vid(m)].pos - m.vertex_attrs[verts[2].vid(m)].pos) + // .cross( + // m.vertex_attrs[verts[1].vid(m)].pos - + // m.vertex_attrs[verts[2].vid(m)].pos)) + // .norm() / + // 2.0; + // }; + // auto normal = [](auto& m, auto& verts) { + // return ((m.vertex_attrs[verts[0].vid(m)].pos - m.vertex_attrs[verts[2].vid(m)].pos) + // .cross( + // m.vertex_attrs[verts[1].vid(m)].pos - + // m.vertex_attrs[verts[2].vid(m)].pos)) + // .normalized(); + // }; + // auto w0 = 0.0; + // Eigen::Vector3d n0(0.0, 0.0, 0.0); + // for (auto& e : one_ring_tris) { + // auto verts = oriented_tri_vertices(e); + // w0 += area(*this, verts); + // n0 += area(*this, verts) * normal(*this, verts); + // } + // n0 /= w0; + // if (n0.norm() < 1e-10) return vertex_attrs[t.vid(*this)].pos; + // n0 = n0.normalized(); + // after_smooth += n0 * n0.transpose() * (vertex_attrs[t.vid(*this)].pos - after_smooth); + // assert(check_mesh_connectivity_validity()); return after_smooth; } -bool UniformRemeshing::uniform_remeshing(double L, int iterations) +bool UniformRemeshing::uniform_remeshing(int iterations, bool debug_output) { - int cnt = 0; - wmtk::logger().info("target len is: {}", L); + if (debug_output) { + update_qualities(); + write_vtu(fmt::format("debug_{}", debug_print_counter++)); + } igl::Timer timer; - while (cnt < iterations) { - cnt++; - wmtk::logger().info("??? Pass ??? {}", cnt); + for (int cnt = 0; cnt < iterations; ++cnt) { + wmtk::logger().info("===== Pass {} of {} =====", cnt, iterations); // split timer.start(); - split_remeshing(L); + split_remeshing(); wmtk::logger().info("--------split time-------: {} ms", timer.getElapsedTimeInMilliSec()); + if (debug_output) { + update_qualities(); + write_vtu(fmt::format("debug_{}", debug_print_counter++)); + } // collpase timer.start(); - collapse_remeshing(L); + collapse_remeshing(); wmtk::logger().info( "--------collapse time-------: {} ms", timer.getElapsedTimeInMilliSec()); + if (debug_output) { + update_qualities(); + write_vtu(fmt::format("debug_{}", debug_print_counter++)); + } // swap edges timer.start(); swap_remeshing(); wmtk::logger().info("--------swap time-------: {} ms", timer.getElapsedTimeInMilliSec()); + if (debug_output) { + update_qualities(); + write_vtu(fmt::format("debug_{}", debug_print_counter++)); + } // smoothing timer.start(); smooth_all_vertices(); wmtk::logger().info("--------smooth time-------: {} ms", timer.getElapsedTimeInMilliSec()); + if (debug_output) { + update_qualities(); + write_vtu(fmt::format("debug_{}", debug_print_counter++)); + } partition_mesh_morton(); } + update_qualities(); wmtk::logger().info("finished {} remeshing iterations", iterations); wmtk::logger().info("+++++++++finished+++++++++"); return true; @@ -707,4 +1213,403 @@ bool UniformRemeshing::write_triangle_mesh(std::string path) assert(manifold); return manifold; -} \ No newline at end of file +} + +double UniformRemeshing::shape_quality(const Vector3d& p0, const Vector3d& p1, const Vector3d& p2) + const +{ + const Eigen::Vector3d a = (p1 - p0); + const Eigen::Vector3d b = (p2 - p1); + const Eigen::Vector3d c = (p0 - p2); + + const double sq_length_sum = a.squaredNorm() + b.squaredNorm() + c.squaredNorm(); + const double area = a.cross(b).norm(); + + const double prefactor = 2 * std::sqrt(3); + + const double quality = (prefactor * area) / sq_length_sum; + + return quality; +} + +double UniformRemeshing::get_quality(const std::array& vs) const +{ + return shape_quality(vertex_attrs[vs[0]].pos, vertex_attrs[vs[1]].pos, vertex_attrs[vs[2]].pos); +} +double UniformRemeshing::get_quality(const Tuple& t) const +{ + const auto vs = oriented_tri_vids(t); + return get_quality(vs); +} + +void UniformRemeshing::update_qualities() +{ + double min_q = std::numeric_limits::max(); + double avg_q = 0; + const auto faces = get_faces(); + for (const Tuple& t : faces) { + const double q = get_quality(t); + face_attrs[t.fid(*this)].quality = q; + min_q = std::min(min_q, q); + avg_q += q; + } + avg_q /= faces.size(); + logger().info("Quality min = {:.2}, avg = {:.2}", min_q, avg_q); +} + +void UniformRemeshing::set_feature_vertices(const std::vector& feature_vertices) +{ + for (size_t i = 0; i < vertex_attrs.size(); ++i) { + vertex_attrs[i].is_feature = false; + } + + for (size_t v : feature_vertices) { + if (v < vertex_attrs.size()) { + vertex_attrs[v].is_feature = true; + } else { + wmtk::logger().warn("set_feature_vertices: index {} out of range", v); + } + } + + wmtk::logger().info("Marked {} feature vertices", feature_vertices.size()); +} + +bool UniformRemeshing::is_feature_vertex(size_t vid) const +{ + return vid < vertex_attrs.size() && vertex_attrs[vid].is_feature; +} + +bool UniformRemeshing::is_feature_edge(const Tuple& t) const +{ + return edge_attrs[t.eid(*this)].is_feature > 0; +} + +void UniformRemeshing::set_feature_edges( + const std::vector, int>>& feature_edges) +{ + m_input_feature_edges = feature_edges; +} + +void UniformRemeshing::set_patch_ids(const std::vector& patch_ids) +{ + if (patch_ids.size() != tri_capacity()) { + log_and_throw_error( + "Patch ID vector has not the right size. Patch IDs {}, #faces {}", + patch_ids.size(), + tri_capacity()); + } + + for (size_t i = 0; i < patch_ids.size(); ++i) { + face_attrs[i].patch_id = patch_ids[i]; + } +} + +void UniformRemeshing::initialize_feature_edges() +{ + logger().info("Init feature edges"); + std::map feature_edges; + + for (const auto& item : m_input_feature_edges) { + const auto& ab = item.first; + const size_t v0 = ab[0]; + const size_t v1 = ab[1]; + if (v0 == v1) { + continue; + } + feature_edges[simplex::Edge(v0, v1)] = item.second; + + if (v0 >= vertex_attrs.size()) { + log_and_throw_error( + "Invalid feature edge ({},{}), vertex id {} is invalid", + v0, + v1, + v0); + } + if (v1 >= vertex_attrs.size()) { + log_and_throw_error( + "Invalid feature edge ({},{}), vertex id {} is invalid", + v0, + v1, + v1); + } + + vertex_attrs[v0].is_feature = true; + vertex_attrs[v1].is_feature = true; + + const size_t eid = eid_from_vids(v0, v1); + edge_attrs[eid].is_feature = item.second; + } + + // size_t found = 0; + // for (const Tuple& e : get_edges()) { + // const size_t v0 = e.vid(*this); + // const size_t v1 = e.switch_vertex(*this).vid(*this); + // + // const simplex::Edge s(v0, v1); + // if (feature_edges.count(s) > 0) { + // const size_t eid = e.eid(*this); + // + // int seg_id = feature_edges[s]; + // // for (const auto& item : m_input_feature_edges) { + // // const auto& ab = item.first; + // // const int id = item.second; + // // if ((ab[0] == v0 && ab[1] == v1) || (ab[0] == v1 && ab[1] == v0)) { + // // seg_id = id; + // // break; + // // } + // // } + // + // edge_attrs[eid].is_feature = seg_id; + // ++found; + // } + //} + + wmtk::logger().info("initialize_feature_edges: marked {} feature edges", feature_edges.size()); + + // if (found != feature_edges.size()) { + // log_and_throw_error( + // "initialize_feature_edges: {} requested feature edges were not found in the mesh", + // feature_edges.size() - found); + // } +} + + +void UniformRemeshing::set_target_edge_length(const double tal) +{ + for (const Tuple& t : get_vertices()) { + vertex_attrs[t.vid(*this)].tal = tal; + } +} + +void UniformRemeshing::set_per_patch_target_edge_length(const double factor) +{ + assert(factor > 0); + + // compute average edge length per patch + std::unordered_map patch_edge_length; + std::unordered_map> patch_edge_lengths; + std::unordered_map patch_face_count; + + for (const Tuple& t : get_faces()) { + const size_t pid = face_attrs[t.fid(*this)].patch_id; + const auto vs = oriented_tri_vids(t); + const Vector3d p0 = vertex_attrs[vs[0]].pos; + const Vector3d p1 = vertex_attrs[vs[1]].pos; + const Vector3d p2 = vertex_attrs[vs[2]].pos; + const double l0 = (p1 - p0).norm(); + const double l1 = (p2 - p1).norm(); + const double l2 = (p0 - p2).norm(); + + if (patch_face_count.count(pid) == 0) { + patch_edge_length[pid] = l0 + l1 + l2; + patch_face_count[pid] = 1; + } else { + patch_edge_length[pid] += l0 + l1 + l2; + patch_face_count[pid] += 1; + } + patch_edge_lengths[pid].emplace_back(l0); + patch_edge_lengths[pid].emplace_back(l1); + patch_edge_lengths[pid].emplace_back(l2); + } + + for (auto& [pid, el] : patch_edge_length) { + //el /= (3 * patch_face_count[pid]); + auto& vec = patch_edge_lengths[pid]; + std::sort(vec.begin(), vec.end()); + //el = *(vec.begin() + std::distance(vec.begin(), vec.end()) / 2); + el = vec.back(); + } + + // apply factor + for (auto& [pid, el] : patch_edge_length) { + el *= factor; + logger().info("Patch {}, tal = {}", pid, el); + } + + // assign shortest patch edge length to vertices + for (const Tuple& t : get_vertices()) { + const auto fs = get_one_ring_tris_for_vertex(t); + double min_el = std::numeric_limits::max(); + for (const Tuple& f : fs) { + const size_t pid = face_attrs[f.fid(*this)].patch_id; + const double patch_el = patch_edge_length[pid]; + min_el = std::min(min_el, patch_el); + } + vertex_attrs[t.vid(*this)].tal = min_el; + } + + // smooth edge length between vertices + for (size_t i = 0; i < 5; ++i) { + for (const Tuple& t : get_vertices()) { + auto& tal = vertex_attrs[t.vid(*this)].tal; + for (const Tuple& tt : get_one_ring_edges_for_vertex(t)) { + tal = std::min(tal, 1.5 * vertex_attrs.at(tt.vid(*this)).tal); + } + } + } +} + +bool UniformRemeshing::write_feature_vertices_obj(const std::string& path) const +{ + std::ofstream out(path); + if (!out.is_open()) return false; + + std::vector fmap; + fmap.reserve(vertex_attrs.size()); + + for (size_t i = 0; i < vertex_attrs.size(); ++i) { + if (!vertex_attrs[i].is_feature) continue; + const auto& p = vertex_attrs[i].pos; + out << "v " << p.x() << " " << p.y() << " " << p.z() << "\n"; + fmap.push_back(i); + } + + for (size_t k = 0; k < fmap.size(); ++k) { + out << "p " << (k + 1) << "\n"; + } + + return true; +} + +void UniformRemeshing::write_vtu(const std::string& path) const +{ + const std::string out_path = path + ".vtu"; + logger().info("Write {}", out_path); + const auto& vs = get_vertices(); + const auto& faces = get_faces(); + + std::vector> edges; + std::vector curve_ids; + for (const Tuple& e : get_edges()) { + auto eid = e.eid(*this); + if (edge_attrs[eid].is_feature) { + size_t v0 = e.vid(*this); + size_t v1 = e.switch_vertex(*this).vid(*this); + edges.emplace_back(std::array{{v0, v1}}); + curve_ids.emplace_back(edge_attrs[eid].is_feature); + } + } + + Eigen::MatrixXd V(vert_capacity(), 3); + Eigen::MatrixXi F(faces.size(), 3); + Eigen::MatrixXi E(edges.size(), 2); + + V.setZero(); + F.setZero(); + E.setZero(); + + Eigen::VectorXd v_is_feature(vert_capacity()); + Eigen::VectorXd v_corner_id(vert_capacity()); + Eigen::VectorXd v_tal(vert_capacity()); + Eigen::VectorXd v_valence(vert_capacity()); + Eigen::VectorXd f_pid(faces.size()); + Eigen::VectorXd f_quality(faces.size()); + Eigen::VectorXd c_id(curve_ids.size()); + + for (size_t i = 0; i < curve_ids.size(); ++i) { + c_id(i) = curve_ids[i] - 1; + } + + int index = 0; + for (const Tuple& t : faces) { + const auto& vs = oriented_tri_vids(t); + for (int j = 0; j < 3; j++) { + F(index, j) = vs[j]; + } + f_pid[index] = face_attrs[t.fid(*this)].patch_id; + f_quality[index] = face_attrs[t.fid(*this)].quality; + ++index; + } + + for (size_t i = 0; i < edges.size(); ++i) { + for (size_t j = 0; j < 2; ++j) { + E(i, j) = edges[i][j]; + } + } + + for (const Tuple& v : vs) { + const size_t vid = v.vid(*this); + V.row(vid) = vertex_attrs[vid].pos; + v_is_feature[vid] = vertex_attrs[vid].is_feature; + v_corner_id[vid] = vertex_attrs[vid].corner_id; + v_tal[vid] = vertex_attrs[vid].tal; + v_valence[vid] = get_valence_for_vertex(vid); + } + + std::shared_ptr writer; + writer = std::make_shared(); + + writer->add_field("is_feature", v_is_feature); + writer->add_field("corner_id", v_corner_id); + writer->add_field("target_edge_length", v_tal); + writer->add_field("valence", v_valence); + writer->add_cell_field("patch_id", f_pid); + writer->add_cell_field("quality", f_quality); + writer->write_mesh(path + ".vtu", V, F); + + // Update report file (if present) without having to pass JSON around. + // Contract: main writes report to `${output}_report.json` where `output` is + // the string passed to `write_vtu(output)`. + const std::string report_file = path + "_report.json"; + + nlohmann::json report = nlohmann::json::object(); + if (std::filesystem::exists(report_file)) { + std::ifstream fin(report_file); + if (fin) { + fin >> report; + } + } + if (!report.is_object()) { + report = nlohmann::json::object(); + } + if (!report.contains("after") || !report["after"].is_object()) { + report["after"] = nlohmann::json::object(); + } + + // min/max internal angle + { + Eigen::VectorXd angles; + igl::internal_angles(V, F, angles); + const auto min_angle = angles.minCoeff(); + const auto max_angle = angles.maxCoeff(); + report["after"]["min_angle"] = min_angle; + report["after"]["max_angle"] = max_angle; + const auto avg_angle = angles.mean(); + report["after"]["avg_angle"] = avg_angle; + } + + // min/max doublearea + { + Eigen::VectorXd double_areas; + igl::doublearea(V, F, double_areas); + const auto min_da = double_areas.minCoeff(); + const auto max_da = double_areas.maxCoeff(); + report["after"]["min_da"] = min_da; + report["after"]["max_da"] = max_da; + // avg double area + const auto avg_da = double_areas.mean(); + report["after"]["avg_da"] = avg_da; + } + + // Persist updated report + { + std::ofstream fout(report_file); + fout << std::setw(4) << report; + } + + // feature edges + { + const auto edge_out_path = path + "_edges.vtu"; + std::shared_ptr edge_writer; + edge_writer = std::make_shared(); + edge_writer->add_field("is_feature", v_is_feature); + edge_writer->add_field("corner_id", v_corner_id); + edge_writer->add_field("target_edge_length", v_tal); + edge_writer->add_field("valence", v_valence); + edge_writer->add_cell_field("curve_id", c_id); + + logger().info("Write {}", edge_out_path); + edge_writer->write_mesh(edge_out_path, V, E); + } +} +} // namespace app::remeshing \ No newline at end of file diff --git a/app/remeshing/src/remeshing/UniformRemeshing.h b/app/remeshing/src/remeshing/UniformRemeshing.h index f9a391e9b1..efefd7d2a7 100644 --- a/app/remeshing/src/remeshing/UniformRemeshing.h +++ b/app/remeshing/src/remeshing/UniformRemeshing.h @@ -1,12 +1,16 @@ #pragma once #include #include +#include #include -#include "wmtk/AttributeCollection.hpp" +#include // clang-format off #include #include +#include +#include +#include #include #include #include @@ -21,6 +25,7 @@ #include #include #include +#include namespace app::remeshing { struct VertexAttributes @@ -28,18 +33,38 @@ struct VertexAttributes Eigen::Vector3d pos; // TODO: in fact, partition id should not be vertex attribute, it is a fixed marker to distinguish tuple/operations. size_t partition_id; - bool freeze = false; + int corner_id = -1; + bool is_feature = false; + double tal = -1; // target edge length +}; + +struct EdgeAttributes +{ + int is_feature = 0; +}; + +struct FaceAttributes +{ + size_t patch_id = -1; + double quality = -1; }; class UniformRemeshing : public wmtk::TriMesh { public: wmtk::SampleEnvelope m_envelope; + wmtk::SampleEnvelope m_feature_envelope; // envelope for feature edges bool m_has_envelope = false; using VertAttCol = wmtk::AttributeCollection; VertAttCol vertex_attrs; + using EdgeAttCol = wmtk::AttributeCollection; + EdgeAttCol edge_attrs; + + using FaceAttCol = wmtk::AttributeCollection; + FaceAttCol face_attrs; + int retry_limit = 10; UniformRemeshing( std::vector _m_vertex_positions, @@ -51,19 +76,78 @@ class UniformRemeshing : public wmtk::TriMesh void create_mesh( size_t n_vertices, const std::vector>& tris, - const std::vector& frozen_verts = std::vector(), + const std::vector>& frozen_verts = {}, bool m_freeze = true, double eps = 0); - struct PositionInfoCache + void initialize_feature_edges(); + + void set_target_edge_length(const double tal); + /** + * @brief Set the target edge length to the per-patch average. + * + * Smooth the edge length in between vertices to have a difference of at max 1.5. + */ + void set_per_patch_target_edge_length(const double factor); + + struct SplitInfoCache + { + // incident vertices + size_t v0 = size_t(-1); + size_t v1 = size_t(-1); + wmtk::Vector3d v0p; + wmtk::Vector3d v1p; + double v0_tal; // v0 target edge length + double v1_tal; // v1 target edge length + + int partition_id; + + int is_feature_edge = 0; + + std::map edge_attrs; + std::unordered_map face_attrs; + }; + tbb::enumerable_thread_specific split_info_cache; + + struct SwapInfoCache { + size_t v0 = size_t(-1); + size_t v1 = size_t(-1); + + // opposite vertices + size_t v2 = size_t(-1); + size_t v3 = size_t(-1); + wmtk::Vector3d v0p; + wmtk::Vector3d v1p; + + int is_feature_edge = 0; + std::map ring_edge_attrs; + }; + tbb::enumerable_thread_specific swap_info_cache; + + + struct CollapseInfoCache + { + Eigen::Vector3d v0p; Eigen::Vector3d v1p; - Eigen::Vector3d v2p; int partition_id; + double v0_tal; // v0 target edge length + double v1_tal; // v1 target edge length + int v0_corner_id; + int v1_corner_id; + + size_t v0 = size_t(-1); + size_t v1 = size_t(-1); + int is_feature_edge = 0; + std::vector> + link_vertices; // store the edge-link vertex v2 and the feature ID of (v0,v2) or (v1,v2) }; - tbb::enumerable_thread_specific position_cache; + tbb::enumerable_thread_specific collapse_info_cache; + + bool cache_edge_positions(const Tuple& t); - void cache_edge_positions(const Tuple& t); + std::vector> get_edges_by_condition( + std::function cond) const; bool invariants(const std::vector& new_tris) override; @@ -101,9 +185,10 @@ class UniformRemeshing : public wmtk::TriMesh bool smooth_before(const Tuple& t) override; bool smooth_after(const Tuple& t) override; - double compute_edge_cost_collapse(const TriMesh::Tuple& t, double L) const; - double compute_edge_cost_split(const TriMesh::Tuple& t, double L) const; - double compute_vertex_valence(const TriMesh::Tuple& t) const; + double compute_edge_cost_collapse(const Tuple& t) const; + double compute_edge_cost_split(const Tuple& t) const; + double compute_vertex_valence(const Tuple& t) const; + double compute_swap_energy(const Tuple& t) const; /** * @brief Report statistics. * @@ -116,11 +201,44 @@ class UniformRemeshing : public wmtk::TriMesh * min valence */ std::vector average_len_valen(); - bool split_remeshing(double L); - bool collapse_remeshing(double L); + bool split_remeshing(); + bool collapse_remeshing(); bool swap_remeshing(); - bool uniform_remeshing(double L, int interations); + bool uniform_remeshing(int interations, bool debug_output = false); bool write_triangle_mesh(std::string path); + + /** + * @brief Compute the shape quality. + * + * The shape quality is the area divided by the sum of the squared edge lengths. + * The quality is normalized, i.e., the best quality is 1 and the worst is 0. + */ + double shape_quality( + const wmtk::Vector3d& p0, + const wmtk::Vector3d& p1, + const wmtk::Vector3d& p2) const; + + /** + * @brief Get the quality of a single triangle. + */ + double get_quality(const std::array& vs) const; + double get_quality(const Tuple& t) const; + /** + * @brief Update quality of all triangles. + */ + void update_qualities(); + + void set_feature_vertices(const std::vector& feature_vertices); + void set_feature_edges(const std::vector, int>>& feature_edges); + void set_patch_ids(const std::vector& patch_ids); + bool is_feature_vertex(size_t vid) const; + bool is_feature_edge(const Tuple& t) const; + bool write_feature_vertices_obj(const std::string& path) const; + + void write_vtu(const std::string& path) const; + +private: + std::vector, int>> m_input_feature_edges; }; -} // namespace app::remeshing +} // namespace app::remeshing \ No newline at end of file diff --git a/app/remeshing/tests/test_remeshing.cpp b/app/remeshing/tests/test_remeshing.cpp index c51945e928..ba8a1a43d5 100644 --- a/app/remeshing/tests/test_remeshing.cpp +++ b/app/remeshing/tests/test_remeshing.cpp @@ -29,7 +29,7 @@ TEST_CASE("uniform_remeshing", "[test_remeshing][.]") for (int j = 0; j < 3; j++) tri[i][j] = (size_t)F(i, j); } UniformRemeshing m(v); - std::vector modified_v; + std::vector> modified_v; m.create_mesh(V.rows(), tri, modified_v, 0); REQUIRE(m.check_mesh_connectivity_validity()); REQUIRE(m.uniform_remeshing(0.01, 5)); @@ -45,10 +45,11 @@ TEST_CASE("split_each_edge", "[test_remeshing]") UniformRemeshing m(v_positions, 0); std::vector> tris = {{{0, 1, 2}}}; - std::vector modified_v; + std::vector> modified_v; m.create_mesh(3, tris, modified_v, 0); + m.set_target_edge_length(1.4); int target_vertnum = m.vert_capacity() + 3 * m.get_edges().size() + 3 * m.tri_capacity(); - m.split_remeshing(2.7 / 2); + m.split_remeshing(); m.consolidate_mesh(); REQUIRE(target_vertnum == 15); @@ -75,7 +76,7 @@ TEST_CASE("test_swap", "[test_remeshing]") for (int j = 0; j < 3; j++) tri[i][j] = (size_t)F(i, j); } UniformRemeshing m(v); - std::vector modified_v; + std::vector> modified_v; m.create_mesh(V.rows(), tri, modified_v, 0); int v_invariant = m.get_vertices().size(); int e_invariant = m.get_edges().size(); @@ -117,72 +118,73 @@ TEST_CASE("test_split", "[test_remeshing]") for (int j = 0; j < 3; j++) tri[i][j] = (size_t)F(i, j); } UniformRemeshing m(v); - std::vector modified_v; + std::vector> modified_v; m.create_mesh(V.rows(), tri, modified_v, 0); - m.split_remeshing(m.average_len_valen()[0] * 0.5); + m.set_target_edge_length(1); + m.split_remeshing(); REQUIRE(m.check_mesh_connectivity_validity()); } -TEST_CASE("remeshing_hanging", "[test_remeshing]") -{ - const std::string root(WMTK_DATA_DIR); - const std::string path = root + "/100071_sf.obj"; - std::string output = "100071_out.obj"; - double env_rel = 1e-3; - double len_rel = 5; - int thread = 0; +// TEST_CASE("remeshing_hanging", "[test_remeshing]") +// { +// const std::string root(WMTK_DATA_DIR); +// const std::string path = root + "/100071_sf.obj"; +// std::string output = "100071_out.obj"; +// double env_rel = 1e-3; +// double len_rel = 5; +// int thread = 0; - wmtk::logger().info("remeshing on {}", path); - Eigen::MatrixXd V; - Eigen::MatrixXi F; - bool ok = igl::read_triangle_mesh(path, V, F); - Eigen::VectorXi SVI, SVJ; - Eigen::MatrixXd temp_V = V; // for STL file - igl::remove_duplicate_vertices(temp_V, 0, V, SVI, SVJ); - for (int i = 0; i < F.rows(); i++) - for (int j : {0, 1, 2}) F(i, j) = SVJ[F(i, j)]; - wmtk::logger().info("Before_vertices#: {} \n Before_tris#: {}", V.rows(), F.rows()); +// wmtk::logger().info("remeshing on {}", path); +// Eigen::MatrixXd V; +// Eigen::MatrixXi F; +// bool ok = igl::read_triangle_mesh(path, V, F); +// Eigen::VectorXi SVI, SVJ; +// Eigen::MatrixXd temp_V = V; // for STL file +// igl::remove_duplicate_vertices(temp_V, 0, V, SVI, SVJ); +// for (int i = 0; i < F.rows(); i++) +// for (int j : {0, 1, 2}) F(i, j) = SVJ[F(i, j)]; +// wmtk::logger().info("Before_vertices#: {} \n Before_tris#: {}", V.rows(), F.rows()); - std::vector v(V.rows()); - std::vector> tri(F.rows()); - for (int i = 0; i < V.rows(); i++) { - v[i] = V.row(i); - } - for (int i = 0; i < F.rows(); i++) { - for (int j = 0; j < 3; j++) tri[i][j] = (size_t)F(i, j); - } +// std::vector v(V.rows()); +// std::vector> tri(F.rows()); +// for (int i = 0; i < V.rows(); i++) { +// v[i] = V.row(i); +// } +// for (int i = 0; i < F.rows(); i++) { +// for (int j = 0; j < 3; j++) tri[i][j] = (size_t)F(i, j); +// } - const Eigen::MatrixXd box_min = V.colwise().minCoeff(); - const Eigen::MatrixXd box_max = V.colwise().maxCoeff(); - const double diag = (box_max - box_min).norm(); - const double envelope_size = env_rel * diag; - Eigen::VectorXi dummy; - std::vector modified_v; - if (!igl::is_edge_manifold(F) || !igl::is_vertex_manifold(F, dummy)) { - auto v1 = v; - auto tri1 = tri; - wmtk::separate_to_manifold(v1, tri1, v, tri, modified_v); - } +// const Eigen::MatrixXd box_min = V.colwise().minCoeff(); +// const Eigen::MatrixXd box_max = V.colwise().maxCoeff(); +// const double diag = (box_max - box_min).norm(); +// const double envelope_size = env_rel * diag; +// Eigen::VectorXi dummy; +// std::vector> modified_v; +// if (!igl::is_edge_manifold(F) || !igl::is_vertex_manifold(F, dummy)) { +// auto v1 = v; +// auto tri1 = tri; +// wmtk::separate_to_manifold(v1, tri1, v, tri, modified_v); +// } - UniformRemeshing m(v, thread); - m.create_mesh(v.size(), tri, modified_v, envelope_size); - REQUIRE(m.check_edge_manifold()); - m.get_vertices(); - std::vector properties = m.average_len_valen(); - wmtk::logger().info( - "edgelen: avg max min valence:avg max min before remesh is: {}", - properties); - wmtk::logger().info("target edge length is: {}", properties[0] * len_rel); - m.uniform_remeshing(properties[0] * len_rel, 2); - m.consolidate_mesh(); - m.write_triangle_mesh(output); - wmtk::logger().info( - "After_vertices#: {} \n\t After_tris#: {}", - m.vert_capacity(), - m.tri_capacity()); - REQUIRE(m.check_edge_manifold()); -} +// UniformRemeshing m(v, thread); +// m.create_mesh(v.size(), tri, modified_v, envelope_size); +// REQUIRE(m.check_edge_manifold()); +// m.get_vertices(); +// std::vector properties = m.average_len_valen(); +// wmtk::logger().info( +// "edgelen: avg max min valence:avg max min before remesh is: {}", +// properties); +// wmtk::logger().info("target edge length is: {}", properties[0] * len_rel); +// m.uniform_remeshing(properties[0] * len_rel, 2); +// m.consolidate_mesh(); +// m.write_triangle_mesh(output); +// wmtk::logger().info( +// "After_vertices#: {} \n\t After_tris#: {}", +// m.vert_capacity(), +// m.tri_capacity()); +// REQUIRE(m.check_edge_manifold()); +// } std::function&)> is_inverted = [](auto& tri) { Eigen::Vector2d a, b, c; @@ -220,7 +222,7 @@ TEST_CASE("operation orient", "[test_remeshing]") for (int j = 0; j < 3; j++) tri[i][j] = (size_t)F(i, j); } UniformRemeshing m(v); - std::vector modified_v; + std::vector> modified_v; m.create_mesh(V.rows(), tri, modified_v, 1); auto fs = m.get_faces(); for (auto f : fs) { @@ -234,7 +236,8 @@ TEST_CASE("operation orient", "[test_remeshing]") REQUIRE(!is_inverted(tri)); } - m.split_remeshing(0.8); + m.set_target_edge_length(1); + m.split_remeshing(); fs = m.get_faces(); for (auto f : fs) { std::array tri = { @@ -247,7 +250,7 @@ TEST_CASE("operation orient", "[test_remeshing]") REQUIRE(!is_inverted(tri)); } - m.collapse_remeshing(1.5); + m.collapse_remeshing(); fs = m.get_faces(); for (auto f : fs) { std::array tri = { diff --git a/src/wmtk/TriMesh.cpp b/src/wmtk/TriMesh.cpp index 6da0732872..55e3800333 100644 --- a/src/wmtk/TriMesh.cpp +++ b/src/wmtk/TriMesh.cpp @@ -352,7 +352,7 @@ bool TriMesh::split_edge(const Tuple& t, std::vector& new_tris) // record the vids that will be modified for roll backs on failure std::vector> old_vertices(2); - std::vector> old_tris(1); + std::vector> old_tris; old_vertices[0] = std::make_pair(vid1, m_vertex_connectivity[vid1]); old_vertices[1] = std::make_pair(vid2, m_vertex_connectivity[vid2]); for (size_t i = 0; i < vid3s.size(); ++i) { @@ -1396,7 +1396,7 @@ TriMesh::Tuple TriMesh::tuple_from_edge(size_t vid1, size_t vid2, size_t fid) co return Tuple(vid1, 3 - (a + b), fid, *this); } -TriMesh::Tuple wmtk::TriMesh::tuple_from_vids(size_t vid0, size_t vid1, size_t vid2) const +TriMesh::Tuple TriMesh::tuple_from_vids(size_t vid0, size_t vid1, size_t vid2) const { const auto& vf0 = m_vertex_connectivity[vid0]; const auto& vf1 = m_vertex_connectivity[vid1]; @@ -1436,17 +1436,46 @@ TriMesh::Tuple wmtk::TriMesh::tuple_from_vids(size_t vid0, size_t vid1, size_t v return Tuple(); } -simplex::Vertex wmtk::TriMesh::simplex_from_vertex(const Tuple& t) const +size_t TriMesh::eid_from_vids(size_t vid0, size_t vid1) const +{ + const std::vector& fids = m_vertex_connectivity[vid0].m_conn_tris; + + // find face that contain m_vid and v_opp + for (const size_t f : fids) { + const auto& f_vids = m_tri_connectivity[f].m_indices; + + int local_v0 = -1; + int local_v1 = -1; + for (size_t i = 0; i < f_vids.size(); ++i) { + if (f_vids[i] == vid0) { + local_v0 = i; + } else if (f_vids[i] == vid1) { + local_v1 = i; + } + } + if (local_v1 == -1) { + continue; + } + assert(local_v0 != -1); + size_t local_eid = 3 - (local_v0 + local_v1); + // fids are sorted --> return smallest fid + return 3 * f + local_eid; + } + + log_and_throw_error("Could not find edge id for vertices ({},{})", vid0, vid1); +} + +simplex::Vertex TriMesh::simplex_from_vertex(const Tuple& t) const { return simplex::Vertex(t.vid(*this)); } -simplex::Edge wmtk::TriMesh::simplex_from_edge(const Tuple& t) const +simplex::Edge TriMesh::simplex_from_edge(const Tuple& t) const { return simplex::Edge(t.vid(*this), t.switch_vertex(*this).vid(*this)); } -simplex::Face wmtk::TriMesh::simplex_from_face(const Tuple& t) const +simplex::Face TriMesh::simplex_from_face(const Tuple& t) const { const auto vs = oriented_tri_vids(t.fid(*this)); return simplex::Face(vs[0], vs[1], vs[2]); diff --git a/src/wmtk/TriMesh.h b/src/wmtk/TriMesh.h index c42fe30508..2e2adace58 100644 --- a/src/wmtk/TriMesh.h +++ b/src/wmtk/TriMesh.h @@ -320,6 +320,7 @@ class TriMesh Tuple tuple_from_edge(size_t vid1, size_t vid2, size_t fid) const; Tuple tuple_from_vids(size_t vid0, size_t vid1, size_t vid2) const; + size_t eid_from_vids(size_t vid0, size_t vid1) const; simplex::Vertex simplex_from_vertex(const Tuple& t) const; simplex::Edge simplex_from_edge(const Tuple& t) const; @@ -600,9 +601,13 @@ class TriMesh * @param t tuple pointing to a vertex * @return one-ring tris number */ + size_t get_valence_for_vertex(const size_t vid) const + { + return m_vertex_connectivity[vid].m_conn_tris.size(); + } size_t get_valence_for_vertex(const Tuple& t) const { - return m_vertex_connectivity[t.vid(*this)].m_conn_tris.size(); + return get_valence_for_vertex(t.vid(*this)); } /**