Added rosy field options

Author: rjc8237

include/holonomy/holonomy/field/frame_field.h

@@ -89,6 +89,48 @@ void write_frame_field(
     const Eigen::MatrixXd& kappa,
     const Eigen::MatrixXi& period_jump);
 
+/**
+ * @brief Write a rosy field to file.
+ * 
+ * The format is
+ * ```
+ * <num_faces>
+ * 4
+ * <dx> <dy> <dz>
+ * ...
+ * ```
+ * where d is a representative direction on the face.
+ * 
+ * @param output_filename: file location to serialize the frame field
+ * @param V: mesh vertices
+ * @param F: mesh faces
+ * @param reference_field: per-face reference tangent direction matrix
+ * @param theta: offset angles of a representative cross field direction relative to the reference
+ */
+void write_rosy_field(
+    const std::string& output_filename,
+    const Eigen::MatrixXd& V,
+    const Eigen::MatrixXi& F,
+    const Eigen::MatrixXd& reference_field,
+    const Eigen::VectorXd& theta);
+
+/**
+ * @brief Load a rosy field from file.
+ * 
+ * The format is
+ * ```
+ * <num_faces>
+ * 4
+ * <dx> <dy> <dz>
+ * ...
+ * ```
+ * where d is a representative direction on the face.
+ * 
+ * @param input_filename: file location of the rosy field
+ * @return per-face representative field direction
+ */
+Eigen::MatrixXd load_rosy_field(const std::string& input_filename);
+
 /**
  * @brief Load a frame field from file.
  * 

include/holonomy/holonomy/field/intrinsic_field.h

@@ -42,6 +42,14 @@ class IntrinsicNRosyField
         Eigen::VectorXd& face_angle,
         Eigen::MatrixXd& corner_kappa,
         Eigen::MatrixXi& corner_period_jump) const;
+    void get_field(
+        const Mesh<Scalar>& m,
+        const std::vector<int>& vtx_reindex,
+        const Eigen::MatrixXi& F,
+        Eigen::VectorXi& reference_corner,
+        Eigen::VectorXd& face_angle,
+        Eigen::MatrixXd& corner_kappa,
+        Eigen::MatrixXi& corner_period_jump) const;
 
     void get_fixed_faces(
         const Mesh<Scalar>& m,
@@ -61,6 +69,16 @@ class IntrinsicNRosyField
         const Eigen::VectorXd& face_theta,
         const Eigen::MatrixXd& corner_kappa,
         const Eigen::MatrixXi& corner_period_jump);
+    void set_field(
+        const Mesh<Scalar>& m,
+        const std::vector<int>& vtx_reindex,
+        const Eigen::MatrixXi& F, 
+        const Eigen::VectorXd& face_theta,
+        const Eigen::MatrixXd& corner_kappa,
+        const Eigen::MatrixXi& corner_period_jump);
+    void set_theta(
+        const Mesh<Scalar>& m,
+        const Eigen::VectorXd& face_theta);
 
     Scalar min_angle = 0.;
 

src/holonomy/field/frame_field.cpp

@@ -73,6 +73,31 @@ Eigen::MatrixXd generate_frame_field(
     return igl::rotate_vectors(reference_field, theta, B1, B2);
 }
 
+void write_rosy_field(
+    const std::string& output_filename,
+    const Eigen::MatrixXd& V,
+    const Eigen::MatrixXi& F,
+    const Eigen::MatrixXd& reference_field,
+    const Eigen::VectorXd& theta)
+{
+    Eigen::MatrixXd rosy_field = generate_frame_field(V, F, reference_field, theta);
+
+    std::ofstream field_file(output_filename, std::ios::out | std::ios::trunc);
+    field_file << F.rows() << std::endl;
+    field_file << "4" << std::endl;
+    for (int f = 0; f < F.rows(); ++f)
+    {
+        for (int j : {0 , 1, 2})
+        {
+            field_file << std::fixed << std::setprecision(17) << rosy_field(f, j) << " ";
+        }
+
+        field_file << std::endl;
+    }
+
+    field_file.close();
+}
+
 void write_frame_field(
     const std::string& output_filename,
     const Eigen::MatrixXd& reference_field,
@@ -114,6 +139,48 @@ void write_frame_field(
     field_file.close();
 }
 
+Eigen::MatrixXd load_rosy_field(const std::string& input_filename)
+{
+    // Open file
+    spdlog::debug("opening field at {}", input_filename);
+    std::ifstream input_file(input_filename);
+    if (!input_file) return {};
+
+    // get number of faces
+    std::string line;
+    std::getline(input_file, line);
+    std::istringstream iss(line);
+    int num_faces;
+    iss >> num_faces;
+    spdlog::debug("{} faces", num_faces);
+    std::getline(input_file, line); // skip start line
+
+    // initialize vectors
+    Eigen::MatrixXd frame_field(num_faces, 3);
+
+    // Read file one face at a time
+    int f = 0;
+    while ((f < num_faces) && (std::getline(input_file, line))) {
+        std::istringstream iss(line);
+        for (int i : {0 , 1, 2})
+        {
+            iss >> frame_field(f, i);
+        }
+
+        ++f;
+    }
+
+    if (num_faces != f)
+    {
+        spdlog::error("Number of faces inconsistent with number of lines");
+    }
+
+    // Close file
+    input_file.close();
+    return frame_field;
+}
+
+
 std::tuple<Eigen::MatrixXd, Eigen::VectorXd, Eigen::MatrixXd, Eigen::MatrixXi>
 load_frame_field(const std::string& filename)
 {

src/holonomy/field/intrinsic_field.cpp

@@ -1686,7 +1686,6 @@ std::vector<Scalar> generate_cones_from_rotation_form_FIXME(
 void IntrinsicNRosyField::compute_principal_matchings(const Mesh<Scalar>& m)
 {
     int num_halfedges = m.n_halfedges();
-    VectorX rotation_form(num_halfedges);
     for (int hij = 0; hij < num_halfedges; ++hij) {
         if (hij < m.opp[hij]) continue;
         int hji = m.opp[hij];
@@ -2551,5 +2550,58 @@ void IntrinsicNRosyField::set_field(
     }
 }
 
+void IntrinsicNRosyField::set_theta(
+    const Mesh<Scalar>& m,
+    const Eigen::VectorXd& face_theta)
+{
+    int num_faces = m.n_faces();
+    for (int fijk = 0; fijk < num_faces; ++fijk)
+    {
+        // get reference halfedge
+        int hij = m.h[fijk];
+
+        // only process original faces
+        if (m.type[hij] == 2) continue;
+
+        // record face angles
+        theta[fijk] = Scalar(face_theta[fijk]);
+
+        // skip doubling on closed mesh
+        if (m.type[hij] == 0) continue;
+
+        // invert theta on reflected face
+        int fjik = m.f[m.R[hij]];
+        theta[fjik] = M_PI - theta[fijk];
+    }
+}
+
+
+void IntrinsicNRosyField::set_field(
+    const Mesh<Scalar>& m,
+    const std::vector<int>& vtx_reindex,
+    const Eigen::MatrixXi& F, 
+    const Eigen::VectorXd& face_theta,
+    const Eigen::MatrixXd& corner_kappa,
+    const Eigen::MatrixXi& corner_period_jump)
+{
+    std::vector<int> face_reindex;
+    arange(F.rows(), face_reindex);
+    set_field(m, vtx_reindex, F, face_reindex, face_theta, corner_kappa, corner_period_jump);
+}
+
+void IntrinsicNRosyField::get_field(
+    const Mesh<Scalar>& m,
+    const std::vector<int>& vtx_reindex,
+    const Eigen::MatrixXi& F, 
+    Eigen::VectorXi& reference_corner,
+    Eigen::VectorXd& face_theta,
+    Eigen::MatrixXd& corner_kappa,
+    Eigen::MatrixXi& corner_period_jump) const
+{
+    std::vector<int> face_reindex;
+    arange(F.rows(), face_reindex);
+    get_field(m, vtx_reindex, F, face_reindex, reference_corner, face_theta, corner_kappa, corner_period_jump);
+}
+
 } // namespace Holonomy
 } // namespace Penner

src/holonomy/interface.cpp

@@ -315,6 +315,14 @@ std::tuple<Mesh<Scalar>, std::vector<int>> generate_mesh(
     const std::vector<Scalar>& Th_hat,
     std::vector<int> free_cones)
 {
+    // set all flat cones if no cones given
+    std::vector<Scalar> Th_hat_temp = Th_hat;
+    if (Th_hat.empty())
+    {
+        Th_hat_temp = std::vector<Scalar>(V.rows(), 2 * PI);
+    }
+
+
     // Convert VF mesh to halfedge
     bool fix_boundary = false;
     std::vector<int> vtx_reindex, indep_vtx, dep_vtx, v_rep, bnd_loops;
@@ -323,7 +331,7 @@ std::tuple<Mesh<Scalar>, std::vector<int>> generate_mesh(
         F,
         uv,
         F_uv,
-        Th_hat,
+        Th_hat_temp,
         vtx_reindex,
         indep_vtx,
         dep_vtx,