#include <algorithm>
#include <fstream>
#include <vector>

// https://github.com/salvipeter/dual-contouring/
#include <dc.hh>

// Note:
// This generates multiple vertices.
// Convert the output to STL to fix this problem.

using namespace DualContouring;

double scale = 8.0;

int main(int argc, char **argv) {
  std::vector<Point3D> points;
  std::ifstream f("gps70.obj");
  f.exceptions(std::ios::failbit | std::ios::badbit);
  while (!f.eof()) {
    char c;
    Point3D p;
    f >> c >> p[0] >> p[1] >> p[2];
    points.push_back(p);
    f >> std::ws;
  }

  // Add base
  for (int x = 0; x < 16; ++x)
    for (int z = 0; z < 9; ++z)
      points.emplace_back(x - 1, -1, z - 1);

  QuadMesh mesh;
  for (auto p : points) {
    p += Vector3D(1, 1, 1);
    mesh.addPoint((p + Vector3D(0, 0, 0)) * scale);
    mesh.addPoint((p + Vector3D(1, 0, 0)) * scale);
    mesh.addPoint((p + Vector3D(1, 1, 0)) * scale);
    mesh.addPoint((p + Vector3D(0, 1, 0)) * scale);
    mesh.addPoint((p + Vector3D(0, 0, 1)) * scale);
    mesh.addPoint((p + Vector3D(1, 0, 1)) * scale);
    mesh.addPoint((p + Vector3D(1, 1, 1)) * scale);
    mesh.addPoint((p + Vector3D(0, 1, 1)) * scale);
  }

  auto addFace = [&](size_t i, size_t j) {
    const auto &p = points[i];
    auto q = p;
    q[j/2] += j % 2 == 0 ? -1 : 1;
    for (const auto &r : points)
      if (q.data == r.data)
        return;
    auto add = [&](size_t a, size_t b, size_t c, size_t d) {
      mesh.addQuad(8*i+a+1, 8*i+b+1, 8*i+c+1, 8*i+d+1);
    };
    switch (j) {
    case 0: add(4,7,3,0); break;
    case 1: add(2,6,5,1); break;
    case 2: add(1,5,4,0); break;
    case 3: add(3,7,6,2); break;
    case 4: add(3,2,1,0); break;
    case 5: add(5,6,7,4); break;
    }
  };

  for (size_t i = 0; i < points.size(); ++i) {
    addFace(i, 0);
    addFace(i, 1);
    addFace(i, 2);
    addFace(i, 3);
    addFace(i, 4);
    addFace(i, 5);
  }
  mesh.writeOBJ("gps70-mesh.obj");
}