diff --git a/.github/workflows/cmake-single-platform.yml b/.github/workflows/cmake-single-platform.yml
new file mode 100644
index 0000000..bb3be7f
--- /dev/null
+++ b/.github/workflows/cmake-single-platform.yml
@@ -0,0 +1,42 @@
+# This starter workflow is for a CMake project running on a single platform. There is a different starter workflow if you need cross-platform coverage.
+# See: https://github.com/actions/starter-workflows/blob/main/ci/cmake-multi-platform.yml
+name: CMake on a single platform
+
+on:
+  push:
+    branches: [ "main" ]
+  pull_request:
+    branches: [ "main" ]
+
+env:
+  # Customize the CMake build type here (Release, Debug, RelWithDebInfo, etc.)
+  BUILD_TYPE: Release
+
+jobs:
+  build:
+    # The CMake configure and build commands are platform agnostic and should work equally well on Windows or Mac.
+    # You can convert this to a matrix build if you need cross-platform coverage.
+    # See: https://docs.github.com/en/free-pro-team@latest/actions/learn-github-actions/managing-complex-workflows#using-a-build-matrix
+    runs-on: ubuntu-latest
+
+    steps:
+    - uses: actions/checkout@v4
+
+    - name: Install deps
+      run: sudo apt install -y libopengl-dev libgl1-mesa-dev libxrandr-dev libxi-dev libxinerama-dev libgl-dev libxcursor-dev libmpfr-dev libmpfrc++-dev
+
+    - name: Configure CMake
+      # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using a single-configuration generator such as make.
+      # See https://cmake.org/cmake/help/latest/variable/CMAKE_BUILD_TYPE.html?highlight=cmake_build_type
+      run: cmake -B ${{github.workspace}}/build -DCMAKE_BUILD_TYPE=${{env.BUILD_TYPE}}
+
+    - name: Build
+      # Build your program with the given configuration
+      run: cmake --build ${{github.workspace}}/build --config ${{env.BUILD_TYPE}}
+
+    - name: Test
+      working-directory: ${{github.workspace}}/build
+      # Execute tests defined by the CMake configuration.
+      # See https://cmake.org/cmake/help/latest/manual/ctest.1.html for more detail
+      run: make test
+
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 2ebc2eb..e31ac9e 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,10 +1,11 @@
-cmake_minimum_required(VERSION 3.25)
+cmake_minimum_required(VERSION 3.30)
 project(PennerOptimization)
 
 set(CMAKE_CXX_STANDARD 17)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 
 add_compile_options(-fPIC)
+add_compile_options(-Wno-register)
 set(DIR_EXT "${CMAKE_CURRENT_SOURCE_DIR}/ext")
 
 # Set options if top level
@@ -33,8 +34,8 @@ endif()
 
 list(PREPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake)
 
+#include(eigen)
 include(boost)
-include(eigen)
 
 # Optionally get multiprecision libraries
 if(USE_MULTIPRECISION OR USE_MULTIPRECISION_OVERLAY)
@@ -92,7 +93,7 @@ if(USE_SUITESPARSE)
 endif()
 
 # Get external libraries
-include(boost)
+#include(boost)
 include(conformal_ideal_delaunay)
 include(libigl)
 include(spectra)
diff --git a/cmake/eigen.cmake b/cmake/eigen.cmake
index 9222bf0..6836812 100644
--- a/cmake/eigen.cmake
+++ b/cmake/eigen.cmake
@@ -31,7 +31,7 @@ else()
     )
     FetchContent_GetProperties(eigen)
     if(NOT eigen_POPULATED)
-        FetchContent_Populate(eigen)
+        FetchContent_MakeAvailable(eigen)
     endif()
     set(EIGEN_INCLUDE_DIRS ${eigen_SOURCE_DIR})
 
diff --git a/cmake/geometry-central.cmake b/cmake/geometry-central.cmake
index d1d5f67..c8609fb 100755
--- a/cmake/geometry-central.cmake
+++ b/cmake/geometry-central.cmake
@@ -7,5 +7,6 @@ FetchContent_Declare(
     geometry-central
     SYSTEM
     GIT_REPOSITORY https://github.com/nmwsharp/geometry-central.git
+    GIT_TAG 03f198748d263cea3faa040c8a399d59c38d67bc
 )
 FetchContent_MakeAvailable(geometry-central)
diff --git a/include/feature/feature/interface.h b/include/feature/feature/interface.h
index 1240651..805c9d8 100644
--- a/include/feature/feature/interface.h
+++ b/include/feature/feature/interface.h
@@ -17,15 +17,18 @@ namespace Feature {
  * @param V: mesh vertices
  * @param F: mesh faces
  * @param use_minimal_forest: (optional) if true, use minimal feature forest instead of spanning tree
+ * @param target_angle: (optional) dihedral angle for features
  * @return refined vertex positions
  * @return refined face indices
  * @return refined feature edges
  * @return refined feature forest edges (spanning or minimal)
+ * @return map from refined faces to parent faces
  */
-std::tuple<Eigen::MatrixXd, Eigen::MatrixXi, std::vector<VertexEdge>, std::vector<VertexEdge>> generate_refined_feature_mesh(
+std::tuple<Eigen::MatrixXd, Eigen::MatrixXi, std::vector<VertexEdge>, std::vector<VertexEdge>, std::vector<int>> generate_refined_feature_mesh(
     const Eigen::MatrixXd& V,
     const Eigen::MatrixXi& F,
-    bool use_minimal_forest=false);
+    bool use_minimal_forest=false,
+    Scalar target_angle=60.);
 
 /**
  * @brief Generate a parameterization aligned to features and the given cross field.
diff --git a/include/feature/feature/surgery/refinement.h b/include/feature/feature/surgery/refinement.h
index fb8fba6..3bcb854 100644
--- a/include/feature/feature/surgery/refinement.h
+++ b/include/feature/feature/surgery/refinement.h
@@ -77,6 +77,8 @@ class IntrinsicRefinementMesh
     Mesh<Scalar>& get_mesh() { return m_mesh; }
     const Mesh<Scalar>& get_mesh() const { return m_mesh; }
 
+    const std::vector<int>& get_face_parent() const { return m_face_parent; }
+
 protected:
     int get_new_face();
     int get_new_vertex();
@@ -86,6 +88,7 @@ class IntrinsicRefinementMesh
 
     Mesh<Scalar> m_mesh;
     std::map<int, std::pair<int, int>> m_endpoints;
+    std::vector<int> m_face_parent;
 
 private:
     virtual int refine_single_face(int face_index);
@@ -99,8 +102,9 @@ class IntrinsicRefinementMesh
  * @return refined vertex positions
  * @return refined face indices
  * @return refined feature edges
+ * @return map from refined faces to parent faces
  */
-std::tuple<Eigen::MatrixXd, Eigen::MatrixXi, std::vector<VertexEdge>> refine_corner_feature_faces(
+std::tuple<Eigen::MatrixXd, Eigen::MatrixXi, std::vector<VertexEdge>, std::vector<int>> refine_corner_feature_faces(
     const FeatureFinder& feature_finder);
 
 /**
@@ -113,8 +117,9 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXi, std::vector<VertexEdge>> refine_cor
  * @return refined face indices
  * @return refined feature edges
  * @return refined feature forest edges (spanning or minimal)
+ * @return map from refined faces to parent faces
  */
-std::tuple<Eigen::MatrixXd, Eigen::MatrixXi, std::vector<VertexEdge>, std::vector<VertexEdge>>
+std::tuple<Eigen::MatrixXd, Eigen::MatrixXi, std::vector<VertexEdge>, std::vector<VertexEdge>, std::vector<int>>
 refine_feature_components(const FeatureFinder& feature_finder, bool use_minimal_forest = false);
 
 } // namespace Feature
diff --git a/include/holonomy/holonomy/field/intrinsic_field.h b/include/holonomy/holonomy/field/intrinsic_field.h
index cfad600..160532d 100644
--- a/include/holonomy/holonomy/field/intrinsic_field.h
+++ b/include/holonomy/holonomy/field/intrinsic_field.h
@@ -72,6 +72,7 @@ class IntrinsicNRosyField
     void remove_greedy_cone_pairs(const Mesh<Scalar>& m);
     void fix_cone_pair(const Mesh<Scalar>& m);
     void fix_zero_cones(const Mesh<Scalar>& m);
+    void move_nearby_cones(const Mesh<Scalar>& m);
     void collapse_adjacent_cones(const Mesh<Scalar>& m);
     void collapse_nearby_cones(const Mesh<Scalar>& m);
 
diff --git a/include/util/util/common.h b/include/util/util/common.h
index 3855744..f8f9c2a 100644
--- a/include/util/util/common.h
+++ b/include/util/util/common.h
@@ -83,6 +83,12 @@ inline void arange(size_t n, std::vector<int>& vec)
     vec.resize(n);
     std::iota(vec.begin(), vec.end(), 0);
 }
+inline std::vector<int> arange(size_t n)
+{
+    std::vector<int> vec;
+    arange(n, vec);
+    return vec;
+}
 
 template <typename OldScalar, typename NewScalar>
 Mesh<NewScalar> change_mesh_type(const Mesh<OldScalar>& m)
diff --git a/py/optimize_impl/energies.py b/py/optimize_impl/energies.py
index c5f0673..0313dae 100644
--- a/py/optimize_impl/energies.py
+++ b/py/optimize_impl/energies.py
@@ -250,6 +250,8 @@ def get_vertex_energy(
     uv_embed[:,:2] = uv[:,:2]
 
     # Get energy
+    if (colormap == 'none'):
+        energy = 0.25 * np.ones(len(v))
     if (colormap == 'scale_factors'):
         energy = best_fit_conformal_vf(v, f, uv_embed, fuv)
     if (colormap == 'sym_dirichlet'):
@@ -293,6 +295,8 @@ def get_face_energy(
         energy = arap_vf(v, f, uv_embed, fuv)
     if (colormap == 'quadratic_sym_dirichlet'):
         energy = quadratic_sym_dirichlet_vf(v, f, uv_embed, fuv)
+    if (colormap == 'none'):
+        energy = 0.25 * np.ones(len(f))
 
     # Optionally use face weighting
     if (use_face_weight):
diff --git a/scripts/holonomy_overlay.py b/scripts/holonomy_overlay.py
index 83d760c..91bae2c 100644
--- a/scripts/holonomy_overlay.py
+++ b/scripts/holonomy_overlay.py
@@ -1,4 +1,4 @@
-# Script to generate an overlay from a similarity mesh
+# Script to generate an overlay from a marked mesh
 
 import os, sys
 script_dir = os.path.dirname(__file__)
@@ -27,10 +27,11 @@ def similarity_overlay_one(args, fname):
     # Skip meshes that are already processed
     try:
         uv_mesh_path = os.path.join(output_dir, name + '_refined_with_uv.obj')
-        V, F = igl.read_triangle_mesh(uv_mesh_path)
-        if (len(V) > 0):
-            print("Skipping processed mesh")
-            return
+        if os.path.isfile(uv_mesh_path):
+            V, F = igl.read_triangle_mesh(uv_mesh_path)
+            if (len(V) > 0):
+                print("Skipping processed mesh")
+                return
     except:
         pass
 
diff --git a/scripts/holonomy_render.py b/scripts/holonomy_render.py
index 1f29566..96a2447 100644
--- a/scripts/holonomy_render.py
+++ b/scripts/holonomy_render.py
@@ -42,7 +42,7 @@ def render_uv_one(args, fname):
     os.makedirs(output_dir, exist_ok=True)
 
     # Get logger
-    log_path = os.path.join(output_dir, name+'_similarity_render.log')
+    log_path = os.path.join(output_dir, name+'_holonomy_render.log')
     logger = script_util.get_logger(log_path)
     logger.info("Rendering {}".format(name))
 
@@ -64,11 +64,6 @@ def render_uv_one(args, fname):
         logger.error("Could not load uv coordinates")
         return
 
-    # Cut mesh along along uv lines
-    #v_cut = script_util.cut_mesh(v3d, f, uv, fuv)
-    #v3d = v_cut
-    #f = fuv
-
     # Need to build double mesh when it has boundary
     is_bd = igl.is_border_vertex(v3d_orig, f_orig)
     build_double = (np.sum(is_bd) != 0)
@@ -99,13 +94,10 @@ def render_uv_one(args, fname):
         blue_size = 0
 
     # Get cut to singularity edges (or build for double meshes)
-    use_boundary_edges = True
     if (args["no_cut"]):
         v_cut_to_sin = []
         f_cut_to_sin = []
         logger.info("Skipping cut to singularity")
-    elif (use_boundary_edges):
-        v_cut_to_sin, f_cut_to_sin = script_util.get_boundary_edges(v3d, uv, f, fuv, bd_thick)
     else:
         try:
             cut_to_sin_list_path = os.path.join(uv_dir, m + "_output", name+"_cut_to_sin_list.pickle")
@@ -155,20 +147,11 @@ def render_uv_one(args, fname):
     u, v = render.get_corner_uv(n, h, to, f, fuv, uv)
 
     # Get colormap for the mesh
-    logger.info('Getting {} colormap'.format(args['colormap']))
-    #r = np.full(len(v3d), 0.25)
-    r = np.full(len(v3d), 0.5)
-    #colormap = cm.get_cmap('YlOrBr')
-    #colormap = cm.get_cmap('OrRd')
-    #colormap = cm.get_cmap('PuRd')
+    if args['render_color'] == 'blue':
+        r = 0.25 + 0.0 * np.abs(r)
+    elif args['render_color'] == 'purple':
+        r = 0.5 + 0.0 * np.abs(r)
     colormap = cm.get_cmap('cool')
-    colormap_scale = args['colormap_scale']
-    logger.info('Using colormap scale {}'.format(colormap_scale))
-    if (colormap_scale > 0):
-        norm = colors.CenteredNorm(colormap_scale*0.5, colormap_scale*0.6)
-    else:
-        norm = colors.CenteredNorm(0, 1)
-        r = r * 0
     norm = colors.NoNorm()
 
     # Render image
@@ -218,7 +201,6 @@ def render_uv_one(args, fname):
 def render_uv_many(args):
     script_util.run_many(render_uv_one, args)
 
-
 def add_render_uv_arguments(parser):
     parser.add_argument(
         "-o", "--output_dir",
@@ -240,20 +222,9 @@ def add_render_uv_arguments(parser):
         default=1
     )
     parser.add_argument(
-        "--colormap",
-        help="energy for colormap",
-        default="scale_factors"
-    )
-    parser.add_argument(
-        "--interpolate_from_original",
-        help="use interpolated colormap from original values",
-        action="store_true"
-    )
-    parser.add_argument(
-        "--colormap_scale",
-        help="ratio to scale colormap by",
-        type=float,
-        default=1
+        "--render_color",
+        help="color for the mesh render",
+        default="blue"
     )
     parser.add_argument(
         "--suffix",
diff --git a/scripts/optimization_scripts/pipeline.py b/scripts/optimization_scripts/pipeline.py
index e25f9f8..bdb38ef 100644
--- a/scripts/optimization_scripts/pipeline.py
+++ b/scripts/optimization_scripts/pipeline.py
@@ -1,6 +1,6 @@
 import os, sys
 script_dir = os.path.dirname(__file__)
-module_dir = os.path.join(script_dir, '..', 'py')
+module_dir = os.path.join(script_dir, '..', '..', 'py')
 sys.path.append(module_dir)
 
 import numpy as np
@@ -100,4 +100,4 @@
             args = script_util.overwrite_args(args_default, args_spec)
 
             # Run method
-            method_func(args)
\ No newline at end of file
+            method_func(args)
diff --git a/src/feature/interface.cpp b/src/feature/interface.cpp
index 05d86b2..3acc28f 100644
--- a/src/feature/interface.cpp
+++ b/src/feature/interface.cpp
@@ -191,24 +191,27 @@ std::tuple<std::vector<Mesh<Scalar>>, std::vector<Eigen::VectorXi>> generate_com
 }
 
 
-std::tuple<Eigen::MatrixXd, Eigen::MatrixXi, std::vector<VertexEdge>, std::vector<VertexEdge>> generate_refined_feature_mesh(
+std::tuple<Eigen::MatrixXd, Eigen::MatrixXi, std::vector<VertexEdge>, std::vector<VertexEdge>, std::vector<int>> generate_refined_feature_mesh(
     const Eigen::MatrixXd& V,
     const Eigen::MatrixXi& F,
-    bool use_minimal_forest)
+    bool use_minimal_forest,
+    Scalar target_angle)
 {
     // generate default features
     FeatureFinder feature_finder(V, F);
-    feature_finder.mark_dihedral_angle_features(60.);
+    feature_finder.mark_dihedral_angle_features(target_angle);
     feature_finder.prune_small_components(4);
     feature_finder.prune_small_features(5);
 
     // refine faces to avoid faces with two boundary edges
-    auto [V_ref_f, F_ref_f, feature_edges_f] = refine_corner_feature_faces(feature_finder);
+    auto [V_ref_f, F_ref_f, feature_edges_f, fn_to_f_f] = refine_corner_feature_faces(feature_finder);
 
     // refine edges to avoid feature components without a forest edge
     feature_finder = FeatureFinder(V_ref_f, F_ref_f);
     feature_finder.mark_features(feature_edges_f);
-    return refine_feature_components(feature_finder, use_minimal_forest);
+    auto [V_ref_e, F_ref_e, feature_edges_e, spanning_edges_e, fn_to_f_e] = refine_feature_components(feature_finder, use_minimal_forest);
+    auto fn_to_f = vector_compose(fn_to_f_f, fn_to_f_e);
+    return std::make_tuple(V_ref_e, F_ref_e, feature_edges_e, spanning_edges_e, fn_to_f);
 }
 
 AlignedMetricGenerator::AlignedMetricGenerator(
diff --git a/src/feature/surgery/refinement.cpp b/src/feature/surgery/refinement.cpp
index 7f10fcf..9707325 100644
--- a/src/feature/surgery/refinement.cpp
+++ b/src/feature/surgery/refinement.cpp
@@ -11,6 +11,7 @@ namespace Feature {
 IntrinsicRefinementMesh::IntrinsicRefinementMesh(const Mesh<Scalar>& m)
     : m_mesh(m)
     , m_endpoints({})
+    , m_face_parent(arange(m.n_faces()))
 {}
 
 // Refine a single face, disregarding symmetry structure and independent vertices
@@ -35,8 +36,12 @@ int IntrinsicRefinementMesh::refine_single_face(int face_index)
     // get new vertex, face and hafledge indices
     int vl = get_new_vertex();
     int fijl = face_index;
+    int face_parent = m_face_parent[face_index];
     int fjkl = get_new_face();
     int fkil = get_new_face();
+    m_face_parent[fijl] = face_parent;
+    m_face_parent[fjkl] = face_parent;
+    m_face_parent[fkil] = face_parent;
     auto [hil, hli] = get_new_edge();
     auto [hjl, hlj] = get_new_edge();
     auto [hkl, hlk] = get_new_edge();
@@ -130,13 +135,19 @@ int IntrinsicRefinementMesh::refine_single_halfedge(int halfedge_index)
     int vj = m.to[hij];
     int vk = m.to[hjk];
     int vl = m.to[hil];
+    int fijk_parent = m_face_parent[fijk];
+    int fjil_parent = m_face_parent[fjil];
 
     // get new vertex, face and hafledge indices
     int vm = get_new_vertex();
     int fmjk = fijk;
     int fmki = get_new_face();
+    m_face_parent[fmjk] = fijk_parent;
+    m_face_parent[fmki] = fijk_parent;
     int fmil = fjil;
     int fmlj = get_new_face();
+    m_face_parent[fmil] = fjil_parent;
+    m_face_parent[fmlj] = fjil_parent;
     int him = hij;
     int hmi = hji;
     auto [hmj, hjm] = get_new_edge();
@@ -351,6 +362,7 @@ int IntrinsicRefinementMesh::get_new_face()
 
     // add null entries for face attributes
     m.h.push_back(-1);
+    m_face_parent.push_back(-1);
 
     return face_index;
 }
@@ -449,7 +461,7 @@ int get_refined_index(int vertex_index, const std::vector<int>& vtx_reindex)
     }
 }
 
-std::tuple<Eigen::MatrixXd, Eigen::MatrixXi, std::vector<VertexEdge>> refine_corner_feature_faces(const FeatureFinder& feature_finder)
+std::tuple<Eigen::MatrixXd, Eigen::MatrixXi, std::vector<VertexEdge>, std::vector<int>> refine_corner_feature_faces(const FeatureFinder& feature_finder)
 {
     const Mesh<Scalar>& m = feature_finder.get_mesh();
     int num_faces = m.n_faces();
@@ -505,7 +517,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXi, std::vector<VertexEdge>> refine_cor
     const auto& V = feature_finder.get_vertex_positions();
     auto [V_ref, F_ref] = refinement_mesh.generate_mesh(V, vtx_reindex);
 
-    return std::make_tuple(V_ref, F_ref, feature_edges);
+    return std::make_tuple(V_ref, F_ref, feature_edges, refinement_mesh.get_face_parent());
 }
 
 
@@ -544,7 +556,7 @@ void make_minimal_forest(const Mesh<Scalar>& m, std::vector<bool>& is_spanning_h
     }
 }
 
-std::tuple<Eigen::MatrixXd, Eigen::MatrixXi, std::vector<VertexEdge>, std::vector<VertexEdge>> refine_feature_components(const FeatureFinder& feature_finder, bool use_minimal_forest)
+std::tuple<Eigen::MatrixXd, Eigen::MatrixXi, std::vector<VertexEdge>, std::vector<VertexEdge>, std::vector<int>> refine_feature_components(const FeatureFinder& feature_finder, bool use_minimal_forest)
 {
     // initialize mesh for refinement
     const Mesh<Scalar>& m = feature_finder.get_mesh();
@@ -635,7 +647,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXi, std::vector<VertexEdge>, std::vecto
     const auto& V = feature_finder.get_vertex_positions();
     auto [V_ref, F_ref] = refinement_mesh.generate_mesh(V, vtx_reindex);
     
-    return std::make_tuple(V_ref, F_ref, feature_edges, spanning_edges);
+    return std::make_tuple(V_ref, F_ref, feature_edges, spanning_edges, refinement_mesh.get_face_parent());
 }
 
 
diff --git a/src/holonomy/field/cross_field.cpp b/src/holonomy/field/cross_field.cpp
index 01357e1..ceb0079 100644
--- a/src/holonomy/field/cross_field.cpp
+++ b/src/holonomy/field/cross_field.cpp
@@ -1,5 +1,5 @@
 #include "holonomy/field/cross_field.h"
-
+/*
 #include <directional/TriMesh.h>
 #include <directional/IntrinsicFaceTangentBundle.h>
 #include <directional/CartesianField.h>
@@ -243,4 +243,5 @@ std::array<Eigen::MatrixXd, 4> reduce_curl(
 }
 
 } // namespace Holonomy
-} // namespace Penner
\ No newline at end of file
+} // namespace Penner
+*/
diff --git a/src/holonomy/field/intrinsic_field.cpp b/src/holonomy/field/intrinsic_field.cpp
index 38c7da7..4ccb0bb 100644
--- a/src/holonomy/field/intrinsic_field.cpp
+++ b/src/holonomy/field/intrinsic_field.cpp
@@ -1800,8 +1800,12 @@ void IntrinsicNRosyField::collapse_adjacent_cones(const Mesh<Scalar>& m)
         if ((!is_boundary[wj]) && (cones[vj] == 4)) continue;
 
         // modify cones
-        spdlog::info("collapsing cones {} and {}", cones[vi], cones[vj]);
         Scalar defect = 4 - cones[vi];
+
+        if ((cones[vi] + defect) < 1) continue;
+        if ((cones[vj] - defect) < 1) continue;
+
+        spdlog::info("collapsing adjacent cones {} and {}", cones[vi], cones[vj]);
         set_period_jump(m, hij, period_jump[hij] - defect);
         cones[vi] += defect;
         cones[vj] -= defect;
@@ -1818,6 +1822,125 @@ void IntrinsicNRosyField::collapse_adjacent_cones(const Mesh<Scalar>& m)
     }
 }
 
+void IntrinsicNRosyField::move_nearby_cones(const Mesh<Scalar>& m)
+{
+    // compute initial cones
+    std::vector<int> cones = generate_cones(m);
+
+    // generate list of boundary vertices
+    int num_vertices = cones.size();
+    std::vector<bool> is_boundary = compute_boundary_vertices(m);
+    std::vector<bool> is_cone(m.n_ind_vertices(), false);
+    std::vector<bool> is_one_away(m.n_ind_vertices(), false);
+    std::vector<bool> is_two_away(m.n_ind_vertices(), false);
+
+    auto adjust_adacent = [&](int wi)
+    {
+        int vi = m.v_rep[wi];
+        is_cone[vi] = true;
+
+        // set one ring
+        int hij = m.out[wi];
+        do
+        {
+            int vj = m.v_rep[m.to[hij]];
+            is_one_away[vj] = true;
+            hij = m.n[m.opp[hij]];
+        } while (hij != m.out[wi]);
+
+        // now set two ring
+        do
+        {
+            int hji = m.opp[hij];
+            int hjk = hji;
+            do
+            {
+                int vk = m.v_rep[m.to[hjk]];
+                if ((vk != vi) && (!is_one_away[vk]))
+                {
+                    is_two_away[vk] = true;
+                }
+                hjk = m.n[m.opp[hjk]];
+            } while (hjk != hji);
+
+            hij = m.n[m.opp[hij]];
+        } while (hij != m.out[wi]);
+    };
+
+    // get all interior cones to place
+    std::vector<std::pair<int, int>> cones_to_place = {};
+    for (int wi = 0; wi < m.n_vertices(); ++wi)
+    {
+        int vi = m.v_rep[wi];
+
+        // check for interior cone
+        if ((cones[vi] == 4) || (is_boundary[wi])) continue;
+        cones_to_place.push_back({wi, 1});
+        is_cone[vi] = true;
+    }
+
+    // place cones, ensuring sufficient distance throughout
+    int num_cones = cones_to_place.size();
+    for (int i = 0; i < num_cones; ++i)
+    {
+        auto [wi, cone_index] = cones_to_place[i];
+        int vi = m.v_rep[wi];
+        int defect = (4 - cones[vi]);
+        is_cone[vi] = false;
+
+        // if adjacent, try to move two away 
+        if (is_one_away[vi])
+        {
+            spdlog::info("moving adjacent cone {}", vi);
+            int hij = m.out[wi];
+            do
+            {
+                int wj =  m.to[hij];
+                int vj =  m.v_rep[wj];
+
+                // check for vertex adjacent to only vi
+                if ((!is_cone[vj]) && (!is_one_away[vj]))
+                {
+                    move_cone(m, vi, vj, defect);
+                    spdlog::info("moved cone to {}", vj);
+                    wi = wj;
+                    vi = vj;
+                    break;
+                }
+
+                hij = m.n[m.opp[hij]];
+            } while (hij != m.out[wi]);
+        }
+
+        // if two away, try to move to free spot
+        if (is_two_away[vi])
+        {
+            spdlog::info("moving two away cone {}", vi);
+            int hij = m.out[wi];
+            do
+            {
+                int wj =  m.to[hij];
+                int vj =  m.v_rep[wj];
+
+                // check for vertex adjacent to only vi
+                if ((!is_cone[vj]) && (!is_one_away[vj]) && (!is_two_away[vj]))
+                {
+                    move_cone(m, vi, vj, defect);
+                    spdlog::info("moved cone to {}", vj);
+                    wi = wj;
+                    vi = vj;
+                    break;
+                }
+
+                hij = m.n[m.opp[hij]];
+            } while (hij != m.out[wi]);
+        }
+
+        // mark region around cone
+        spdlog::info("placed cone at {}", vi);
+        adjust_adacent(wi);
+    }
+}
 
 void IntrinsicNRosyField::collapse_nearby_cones(const Mesh<Scalar>& m)
 {
@@ -1879,19 +2002,36 @@ void IntrinsicNRosyField::collapse_nearby_cones(const Mesh<Scalar>& m)
             // modify cones
             if (adjacent_cones.size() > 1)
             {
-                int hij = adjacent_cones.front();
+                // get arbitrary cone to remove
                 int hik = adjacent_cones.back();
-                int vj =  m.v_rep[m.to[hij]];
                 int vk =  m.v_rep[m.to[hik]];
-                spdlog::info("collapsing cones {} and {}", cones[vj], cones[vk]);
+                adjacent_cones.pop_back();
+
+                // get defect for removing vk
                 Scalar defect = 4 - cones[vk];
                 if (is_boundary[m.to[hik]]) defect = 2 - cones[vk];
-                set_period_jump(m, hij, period_jump[hij] - defect);
-                set_period_jump(m, hik, period_jump[hik] + defect);
-                cones[vk] += defect;
-                cones[vj] -= defect;
-                spdlog::info("new cones {} and {}", cones[vj], cones[vk]);
-                collapsed = true;
+                while (true)
+                {
+                    if (adjacent_cones.empty()) break;
+
+                    int hij = adjacent_cones.back();
+                    adjacent_cones.pop_back();
+                    int vj =  m.v_rep[m.to[hij]];
+
+                    // cannot collapse cone to below pi/2
+                    if ((cones[vk] + defect) < 1) continue;
+                    if ((cones[vj] - defect) < 1) continue;
+
+                    // collapse cone
+                    spdlog::info("collapsing nearby cones {} and {}", cones[vj], cones[vk]);
+                    set_period_jump(m, hij, period_jump[hij] - defect);
+                    set_period_jump(m, hik, period_jump[hik] + defect);
+                    cones[vk] += defect;
+                    cones[vj] -= defect;
+                    spdlog::info("new cones {} and {}", cones[vj], cones[vk]);
+                    collapsed = true;
+                    break;
+                }
             }
                 
         }
@@ -1900,6 +2040,10 @@ void IntrinsicNRosyField::collapse_nearby_cones(const Mesh<Scalar>& m)
         std::vector<int> final_cones = generate_cones(m);
         for (int vi = 0; vi < num_vertices; ++vi)
         {
+            if (cones[vi] == 0)
+            {
+                spdlog::error("zero cone at {}", vi);
+            }
             if (cones[vi] != final_cones[vi])
             {
                 spdlog::error("inconsistent cones {} and {}", cones[vi], final_cones[vi]);
diff --git a/src/holonomy/interface.cpp b/src/holonomy/interface.cpp
index 814e58f..b57e938 100644
--- a/src/holonomy/interface.cpp
+++ b/src/holonomy/interface.cpp
@@ -407,8 +407,8 @@ std::tuple<VectorX, std::vector<Scalar>> generate_intrinsic_rotation_form(
         use_discrete_metric);
 
     // compute rotation form
-    VectorX rotation_form = generate_intrinsic_rotation_form(*cone_metric, field_params);
     //VectorX rotation_form = generate_intrinsic_rotation_form(*cone_metric, vtx_reindex, V, field_params);
+    VectorX rotation_form = generate_intrinsic_rotation_form(*cone_metric, field_params);
 
     // generate cones from the rotation form
     bool has_bd = (cone_metric->type[0] != 0);
@@ -631,4 +631,4 @@ void optimize_triangle_quality(MarkedPennerConeMetric& marked_metric, double max
 }
 
 } // namespace Holonomy
-} // namespace Penner
\ No newline at end of file
+} // namespace Penner
diff --git a/src/optimization/parameterization/layout.cpp b/src/optimization/parameterization/layout.cpp
index 4d97bfc..330ee22 100644
--- a/src/optimization/parameterization/layout.cpp
+++ b/src/optimization/parameterization/layout.cpp
@@ -1038,7 +1038,7 @@ get_consistent_layout(
       int current_seg = m_o.first_segment[hij];
       while (m_o.vertex_type[m_o.to[current_seg]] != ORIGINAL_VERTEX)
       {
-        std::vector<Scalar> tmp = m_o.seg_bcs[current_seg];
+        std::array<Scalar, 2> tmp = m_o.seg_bcs[current_seg];
         lambdas.push_back(tmp[1]);
         current_seg = m_o.next_segment(current_seg);
       }
@@ -1048,7 +1048,7 @@ get_consistent_layout(
       int cnt = lambdas.size() - 1;
       while (m_o.vertex_type[m_o.to[current_seg]] != ORIGINAL_VERTEX)
       {
-        std::vector<Scalar> tmp = m_o.seg_bcs[current_seg];
+        std::array<Scalar, 2> tmp = m_o.seg_bcs[current_seg];
         if (abs(lambdas[cnt] - tmp[0]) > 1e-10)
         {
           spdlog::error("alignment problem: {}, {}", lambdas[cnt], tmp[0]);