#include <algorithm>
#include <cmath>
#include <iterator>
#include <vector>

#include "dc.hh"

using namespace DualContouring;

inline double norm(const Vector3D &v) {
  return std::sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
}

inline Vector3D normalize(const Vector3D &v) {
  double d = norm(v);
  return d > 0 ? v / d : v;
}

inline double dot(const Vector3D &u, const Vector3D &v) {
  return u[0] * v[0] + u[1] * v[1] + u[2] * v[2];
}

auto halfspace(const Point3D &point, const Vector3D &normal) {
  return [=](const Point3D &p) {
    return dot(p - point, normal);
  };
}

auto cone(double radius, const Vector3D &direction, double aperture) {
  aperture *= M_PI / 180;
  Point3D apex = direction * (1 + radius / tan(aperture));
  return [=](const Point3D &p) {
    auto u = normalize(p - apex);
    double alpha = std::acos(std::min(1.0, std::abs(dot(u, direction))));
    return std::sin(alpha - aperture) * norm(p - apex);
  };
};

// Elliptic blend
auto blend(const std::vector<std::function<double(const Point3D &)>> &fs,
           double range, int exponent) {
  return [=](const Point3D &p) {
    double result = 1;
    for (const auto &f : fs)
      result -= std::pow(std::max(0.0, 1 - f(p) / range), exponent);
    return result;
  };
}

auto leg(const Vector3D &v, double r, double h) {
  auto c1 = cone(r, v, 12);
  auto c2 = cone(r, v, 45);
  auto p1 = halfspace(v * h, v);
  auto p2 = halfspace({}, v * -1);
  return [=](const Point3D &p) {
    return std::max(c1(p), std::max(c2(p), std::max(p1(p), p2(p))));
  };
}

auto tetrapod(double r, double h) {
  double r3 = std::sqrt(3);
  auto l1 = leg(Vector3D(-1,  1, -1) / r3, r, h);
  auto l2 = leg(Vector3D( 1, -1, -1) / r3, r, h);
  auto l3 = leg(Vector3D( 1,  1,  1) / r3, r, h);
  auto l4 = leg(Vector3D(-1, -1,  1) / r3, r, h);
  return blend({ l1, l2, l3, l4 }, 0.03, 5);
  // return [=](const Point3D &p) {
  //   return std::min(l1(p), std::min(l2(p), std::min(l3(p), l4(p))));
  // };
}

int main(int argc, char **argv) {
  size_t res = 300;
  auto t1 = tetrapod(0.3, 1.07);
  auto t2 = tetrapod(0.25, 1.02);
  auto f = [&](const Point3D &p) {
    return std::max(t1(p), -t2(p)); // for FDM printing
//  return t1(p);                   // for SLA printing
  };
  auto mesh = DualContouring::isosurface(f, 0,
                                         { { { -1.2, -1.2, -1.2 },
                                             {  1.2,  1.2,  1.2 } } },
                                         { res, res, res });
  mesh.writeOBJ("tetrapod.obj");
}