diff options
Diffstat (limited to 'extern/quadriflow/src/post-solver.cpp')
| -rw-r--r-- | extern/quadriflow/src/post-solver.cpp | 427 |
1 files changed, 427 insertions, 0 deletions
diff --git a/extern/quadriflow/src/post-solver.cpp b/extern/quadriflow/src/post-solver.cpp new file mode 100644 index 00000000000..6027ddd2eb7 --- /dev/null +++ b/extern/quadriflow/src/post-solver.cpp @@ -0,0 +1,427 @@ +// +// post-solver.cpp +// parametrize +// +// Created by Jingwei on 2/5/18. +// +#include <algorithm> +#include <boost/program_options.hpp> +#include <cmath> +#include <cstdio> +#include <string> + +#include "ceres/ceres.h" +#include "ceres/rotation.h" + +#include "post-solver.hpp" +#include "serialize.hpp" + +namespace qflow { + +/// Coefficient of area constraint. The magnitude is 1 if area is equal to 0. +const double COEFF_AREA = 1; +/// Coefficient of tangent constraint. The magnitude is 0.03 if the bais is reference_length. +/// This is because current tangent constraint is not very accurate. +/// This optimization conflicts with COEFF_AREA. +const double COEFF_TANGENT = 0.02; +/// Coefficient of normal constraint. The magnitude is the arc angle. +const double COEFF_NORMAL = 1; +/// Coefficient of normal constraint. The magnitude is the arc angle. +const double COEFF_FLOW = 1; +/// Coefficient of orthogonal edge. The magnitude is the arc angle. +const double COEFF_ORTH = 1; +/// Coefficient of edge length. The magnitude is the arc angle. +const double COEFF_LENGTH = 1; +/// Number of iterations of the CGNR solver +const int N_ITER = 100; + +template <typename T, typename T2> +T DotProduct(const T a[3], const T2 b[3]) { + return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; +} + +template <typename T> +T Length2(const T a[3]) { + return DotProduct(a, a); +} + +namespace ceres { +inline double min(const double f, const double g) { return std::min(f, g); } + +template <typename T, int N> +inline Jet<T, N> min(const Jet<T, N>& f, const Jet<T, N>& g) { + if (f.a < g.a) + return f; + else + return g; +} +} // namespace ceres + +bool DEBUG = 0; +struct FaceConstraint { + FaceConstraint(double coeff_area, double coeff_normal, double coeff_flow, double coeff_orth, + double length, Vector3d normal[4], Vector3d Q0[4], Vector3d Q1[4]) + : coeff_area(coeff_area), + coeff_normal(coeff_normal), + coeff_flow(coeff_flow), + coeff_orth(coeff_orth), + area0(length * length), + normal0{ + normal[0], + normal[1], + normal[2], + normal[3], + }, + Q0{Q0[0], Q0[1], Q0[2], Q0[3]}, + Q1{Q1[0], Q1[1], Q1[2], Q1[3]} {} + + template <typename T> + bool operator()(const T* p0, const T* p1, const T* p2, const T* p3, T* r) const { + const T* p[] = {p0, p1, p2, p3}; + r[12] = T(); + for (int k = 0; k < 4; ++k) { + auto pc = p[k]; + auto pa = p[(k + 1) % 4]; + auto pb = p[(k + 3) % 4]; + + T a[3]{pa[0] - pc[0], pa[1] - pc[1], pa[2] - pc[2]}; + T b[3]{pb[0] - pc[0], pb[1] - pc[1], pb[2] - pc[2]}; + + T length_a = ceres::sqrt(Length2(a)); + T length_b = ceres::sqrt(Length2(b)); + T aa[3]{a[0] / length_a, a[1] / length_a, a[2] / length_a}; + T bb[3]{b[0] / length_b, b[1] / length_b, b[2] / length_b}; + r[3 * k + 0] = coeff_orth * DotProduct(aa, bb); + + T degree_edge0 = ceres::abs(DotProduct(aa, &Q0[k][0])); + T degree_edge1 = ceres::abs(DotProduct(aa, &Q1[k][0])); + T degree_edge = ceres::min(degree_edge0, degree_edge1); + r[3 * k + 1] = coeff_flow * degree_edge; + + T normal[3]; + ceres::CrossProduct(a, b, normal); + T area = ceres::sqrt(Length2(normal)); + r[12] += area; + + assert(area != T()); + for (int i = 0; i < 3; ++i) normal[i] /= area; + T degree_normal = DotProduct(normal, &normal0[k][0]) - T(1); + r[3 * k + 2] = coeff_normal * degree_normal * degree_normal; + } + r[12] = coeff_area * (r[12] / (4.0 * area0) - 1.0); + return true; + } + + static ceres::CostFunction* create(double coeff_area, double coeff_normal, double coeff_flow, + double coeff_orth, double length, Vector3d normal[4], + Vector3d Q0[4], Vector3d Q1[4]) { + return new ceres::AutoDiffCostFunction<FaceConstraint, 13, 3, 3, 3, 3>(new FaceConstraint( + coeff_area, coeff_normal, coeff_flow, coeff_orth, length, normal, Q0, Q1)); + } + + double coeff_area; + double coeff_normal; + double coeff_flow; + double coeff_orth; + + double area0; + Vector3d normal0[4]; + Vector3d Q0[4], Q1[4]; +}; + +struct VertexConstraint { + VertexConstraint(double coeff_tangent, Vector3d normal, double bias, double length) + : coeff{coeff_tangent / length * 10}, bias0{bias}, normal0{normal} {} + + template <typename T> + bool operator()(const T* p, T* r) const { + r[0] = coeff * (DotProduct(p, &normal0[0]) - bias0); + return true; + } + + static ceres::CostFunction* create(double coeff_tangent, Vector3d normal, double bias, + double length) { + return new ceres::AutoDiffCostFunction<VertexConstraint, 1, 3>( + new VertexConstraint(coeff_tangent, normal, bias, length)); + } + + double coeff; + double bias0; + Vector3d normal0; +}; + +void solve(std::vector<Vector3d>& O_quad, std::vector<Vector3d>& N_quad, + std::vector<Vector3d>& Q_quad, std::vector<Vector4i>& F_quad, + std::vector<double>& B_quad, MatrixXd& V, MatrixXd& N, MatrixXd& Q, MatrixXd& O, + MatrixXi& F, double reference_length, double coeff_area, double coeff_tangent, + double coeff_normal, double coeff_flow, double coeff_orth) { + printf("Parameter: \n"); + printf(" coeff_area: %.4f\n", coeff_area); + printf(" coeff_tangent: %.4f\n", coeff_tangent); + printf(" coeff_normal: %.4f\n", coeff_normal); + printf(" coeff_flow: %.4f\n", coeff_flow); + printf(" coeff_orth: %.4f\n\n", coeff_orth); + int n_quad = Q_quad.size(); + + ceres::Problem problem; + std::vector<double> solution(n_quad * 3); + for (int vquad = 0; vquad < n_quad; ++vquad) { + solution[3 * vquad + 0] = O_quad[vquad][0]; + solution[3 * vquad + 1] = O_quad[vquad][1]; + solution[3 * vquad + 2] = O_quad[vquad][2]; + } + + // Face constraint (area and normal direction) + for (int fquad = 0; fquad < F_quad.size(); ++fquad) { + auto v = F_quad[fquad]; + Vector3d normal[4], Q0[4], Q1[4]; + for (int k = 0; k < 4; ++k) { + normal[k] = N_quad[v[k]]; + Q0[k] = Q_quad[v[k]]; + Q1[k] = Q0[k].cross(normal[k]).normalized(); + } + ceres::CostFunction* cost_function = FaceConstraint::create( + coeff_area, coeff_normal, coeff_flow, coeff_orth, reference_length, normal, Q0, Q1); + problem.AddResidualBlock(cost_function, nullptr, &solution[3 * v[0]], &solution[3 * v[1]], + &solution[3 * v[2]], &solution[3 * v[3]]); + } + + // Tangent constraint + for (int vquad = 0; vquad < O_quad.size(); ++vquad) { + ceres::CostFunction* cost_function = VertexConstraint::create( + coeff_tangent, N_quad[vquad], B_quad[vquad], reference_length); + problem.AddResidualBlock(cost_function, nullptr, &solution[3 * vquad]); + } + + // Flow constraint + + ceres::Solver::Options options; + options.num_threads = 1; + options.max_num_iterations = N_ITER; + options.initial_trust_region_radius = 1; + options.linear_solver_type = ceres::CGNR; + options.minimizer_progress_to_stdout = true; + ceres::Solver::Summary summary; + ceres::Solve(options, &problem, &summary); + + std::cout << summary.BriefReport() << std::endl; + + for (int vquad = 0; vquad < n_quad; ++vquad) { + O_quad[vquad][0] = solution[3 * vquad + 0]; + O_quad[vquad][1] = solution[3 * vquad + 1]; + O_quad[vquad][2] = solution[3 * vquad + 2]; + } + + return; +} + +void optimize_quad_positions(std::vector<Vector3d>& O_quad, std::vector<Vector3d>& N_quad, + std::vector<Vector3d>& Q_quad, std::vector<Vector4i>& F_quad, + VectorXi& V2E_quad, std::vector<int>& E2E_quad, MatrixXd& V, + MatrixXd& N, MatrixXd& Q, MatrixXd& O, MatrixXi& F, VectorXi& V2E, + VectorXi& E2E, DisajointTree& disajoint_tree, double reference_length, + bool just_serialize) { + printf("Quad mesh info:\n"); + printf("Number of vertices with normals and orientations: %d = %d = %d\n", (int)O_quad.size(), + (int)N_quad.size(), (int)Q_quad.size()); + printf("Number of faces: %d\n", (int)F_quad.size()); + printf("Number of directed edges: %d\n", (int)E2E_quad.size()); + // Information for the original mesh + printf("Triangle mesh info:\n"); + printf( + "Number of vertices with normals, " + "orientations and associated quad positions: " + "%d = %d = %d = %d\n", + (int)V.cols(), (int)N.cols(), (int)Q.cols(), (int)O.cols()); + printf("Number of faces: %d\n", (int)F.cols()); + printf("Number of directed edges: %d\n", (int)E2E.size()); + printf("Reference length: %.2f\n", reference_length); + + int flip_count = 0; + for (int i = 0; i < F_quad.size(); ++i) { + bool flipped = false; + for (int j = 0; j < 4; ++j) { + int v1 = F_quad[i][j]; + int v2 = F_quad[i][(j + 1) % 4]; + int v3 = F_quad[i][(j + 3) % 4]; + + Vector3d face_norm = (O_quad[v2] - O_quad[v1]).cross(O_quad[v3] - O_quad[v1]); + Vector3d vertex_norm = N_quad[v1]; + if (face_norm.dot(vertex_norm) < 0) { + flipped = true; + } + } + if (flipped) { + flip_count++; + } + } + printf("Flipped Quads: %d\n", flip_count); + + int n_quad = O_quad.size(); + int n_trig = O.cols(); + std::vector<double> B_quad(n_quad); // Average bias for quad vertex + std::vector<int> B_weight(n_quad); + + printf("ntrig: %d, disjoint_tree.size: %d\n", n_trig, (int)disajoint_tree.indices.size()); + for (int vtrig = 0; vtrig < n_trig; ++vtrig) { + int vquad = disajoint_tree.Index(vtrig); + double b = N_quad[vquad].dot(O.col(vtrig)); + B_quad[vquad] += b; + B_weight[vquad] += 1; + } + for (int vquad = 0; vquad < n_quad; ++vquad) { + assert(B_weight[vquad]); + B_quad[vquad] /= B_weight[vquad]; + } + + puts("Save parameters to post.bin for optimization"); + FILE* out = fopen("post.bin", "wb"); + assert(out); + Save(out, O_quad); + Save(out, N_quad); + Save(out, Q_quad); + Save(out, F_quad); + Save(out, B_quad); + Save(out, V); + Save(out, N); + Save(out, Q); + Save(out, O); + Save(out, F); + Save(out, reference_length); + fclose(out); + + if (!just_serialize) { + puts("Start post optimization"); + solve(O_quad, N_quad, Q_quad, F_quad, B_quad, V, N, Q, O, F, reference_length, COEFF_AREA, + COEFF_TANGENT, COEFF_NORMAL, COEFF_FLOW, COEFF_ORTH); + } +} + +#ifdef POST_SOLVER + +void SaveObj(const std::string& fname, std::vector<Vector3d> O_quad, + std::vector<Vector4i> F_quad) { + std::ofstream os(fname); + for (int i = 0; i < (int)O_quad.size(); ++i) { + os << "v " << O_quad[i][0] << " " << O_quad[i][1] << " " << O_quad[i][2] << "\n"; + } + for (int i = 0; i < (int)F_quad.size(); ++i) { + os << "f " << F_quad[i][0] + 1 << " " << F_quad[i][1] + 1 << " " << F_quad[i][2] + 1 << " " + << F_quad[i][3] + 1 << "\n"; + } + os.close(); +} + +int main(int argc, char* argv[]) { + double coeff_area; + double coeff_tangent; + double coeff_normal; + double coeff_flow; + double coeff_orth; + + namespace po = boost::program_options; + po::options_description desc("Allowed options"); + desc.add_options() // clang-format off + ("help,h", "produce help message") + ("area,a", po::value<double>(&coeff_area)->default_value(COEFF_AREA), "Set the coefficient of area constraint") + ("tangent,t", po::value<double>(&coeff_tangent)->default_value(COEFF_TANGENT), "Set the coefficient of tangent constraint") + ("normal,n", po::value<double>(&coeff_normal)->default_value(COEFF_NORMAL), "Set the coefficient of normal constraint") + ("flow,f", po::value<double>(&coeff_flow)->default_value(COEFF_FLOW), "Set the coefficient of flow (Q) constraint") + ("orth,o", po::value<double>(&coeff_orth)->default_value(COEFF_ORTH), "Set the coefficient of orthogonal constraint"); + + // clang-format on + po::variables_map vm; + po::store(po::parse_command_line(argc, argv, desc), vm); + po::notify(vm); + if (vm.count("help")) { + std::cout << desc << std::endl; + return 1; + } + + std::vector<Vector3d> O_quad; + std::vector<Vector3d> N_quad; + std::vector<Vector3d> Q_quad; + std::vector<Vector4i> F_quad; + std::vector<double> B_quad; + MatrixXd V; + MatrixXd N; + MatrixXd Q; + MatrixXd O; + MatrixXi F; + double reference_length; + + puts("Read parameters from post.bin"); + FILE* in = fopen("post.bin", "rb"); + assert(in); + Read(in, O_quad); + Read(in, N_quad); + Read(in, Q_quad); + Read(in, F_quad); + Read(in, B_quad); + Read(in, V); + Read(in, N); + Read(in, Q); + Read(in, O); + Read(in, F); + Read(in, reference_length); + fclose(in); + printf("reference_length: %.2f\n", reference_length); + SaveObj("presolver.obj", O_quad, F_quad); + + int n_flip = 0; + double sum_degree = 0; + for (int i = 0; i < F_quad.size(); ++i) { + bool flipped = false; + for (int j = 0; j < 4; ++j) { + int v1 = F_quad[i][j]; + int v2 = F_quad[i][(j + 1) % 4]; + int v3 = F_quad[i][(j + 3) % 4]; + + Vector3d face_norm = + (O_quad[v2] - O_quad[v1]).cross(O_quad[v3] - O_quad[v1]).normalized(); + Vector3d vertex_norm = N_quad[v1]; + if (face_norm.dot(vertex_norm) < 0) { + flipped = true; + } + double degree = std::acos(face_norm.dot(vertex_norm)); + assert(degree >= 0); + // printf("cos theta = %.2f\n", degree); + sum_degree += degree * degree; + } + n_flip += flipped; + } + printf("n_flip: %d\nsum_degree: %.3f\n", n_flip, sum_degree); + + puts("Start post optimization"); + solve(O_quad, N_quad, Q_quad, F_quad, B_quad, V, N, Q, O, F, reference_length, coeff_area, + coeff_tangent, coeff_normal, coeff_flow, coeff_orth); + SaveObj("postsolver.obj", O_quad, F_quad); + + n_flip = 0; + sum_degree = 0; + for (int i = 0; i < F_quad.size(); ++i) { + bool flipped = false; + for (int j = 0; j < 4; ++j) { + int v1 = F_quad[i][j]; + int v2 = F_quad[i][(j + 1) % 4]; + int v3 = F_quad[i][(j + 3) % 4]; + + Vector3d face_norm = + (O_quad[v2] - O_quad[v1]).cross(O_quad[v3] - O_quad[v1]).normalized(); + Vector3d vertex_norm = N_quad[v1]; + if (face_norm.dot(vertex_norm) < 0) { + flipped = true; + } + double degree = std::acos(face_norm.dot(vertex_norm)); + assert(degree >= 0); + sum_degree += degree * degree; + } + n_flip += flipped; + } + printf("n_flip: %d\nsum_degree: %.3f\n", n_flip, sum_degree); + return 0; +} + +#endif + +} // namespace qflow |