quadriflow/src/post-solver.cpp

//
//  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