#include <OpenMesh/Core/Mesh/PolyMesh_ArrayKernelT.hh>
#include <OpenMesh/Core/IO/MeshIO.hh>
#include <OpenMesh/Tools/Subdivider/Uniform/CatmullClarkT.hh>

// Parameters
size_t iterations = 100;
size_t subd_iterations = 5;

using Mesh = OpenMesh::PolyMesh_ArrayKernelT<OpenMesh::DefaultTraits>;
using Vector = OpenMesh::Vec3f;

void subdivide(Mesh &m, int n) {
  OpenMesh::Subdivider::Uniform::CatmullClarkT<Mesh> cc;
  cc.attach(m);
  cc(n);
  cc.detach();
}

int main(int argc, char **argv) {
  if (argc != 2) {
    std::cerr << "Usage: " << argv[0] << " <input-mesh>" << std::endl;
    return 1;
  }

  Mesh original;
  if (!OpenMesh::IO::read_mesh(original, argv[1])) {
    std::cerr << "Cannot read mesh from file: " << argv[1] << std::endl;
    return 2;
  }
  Mesh cage = original;

  for (size_t iter = 1; iter <= iterations; iter++) {
    std::cout << "Iteration " << iter << std::endl;
    Mesh surf = cage;
    subdivide(surf, subd_iterations);
    float maxd = 0;
    for (auto v : cage.vertices()) {
      auto p_orig = original.point(v);
      auto p_cage = cage.point(v);
      auto p_surf = surf.point(v);
      auto d = p_orig - p_surf;
      cage.set_point(v, p_cage + d);
      if (d.norm() > maxd)
        maxd = d.norm();
    }
    std::cout << "Max. change: " << maxd << std::endl;
  }

  OpenMesh::IO::write_mesh(cage, "cage.obj");
  subdivide(cage, subd_iterations);
  OpenMesh::IO::write_mesh(cage, "surface.obj");
}