#include "fiko/dokum/voxelize.hpp"

#include <algorithm>
#include <cmath>
#include <thread>

#include "fiko/geometri/brep.hpp"
#include "fiko/render/mesh_veri.hpp"

namespace fiko::dokum {

static inline Vec3d capraz(Vec3d a, Vec3d b) {
  return {a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x};
}

UcgenCorba corbaBrepten(const geometri::BrepGovde& govde, f64 kaliteOrani) {
  UcgenCorba c;
  if (!govde.gecerliMi()) return c;
  const MeshVeri m = govde.ucgenle({0.6f, 0.6f, 0.6f}, kaliteOrani);
  c.vertler.reserve(m.indeksler.size());
  for (size_t i = 0; i + 2 < m.indeksler.size(); i += 3) {
    auto v = [&](u32 idx) {
      const auto& p = m.vertexler[idx].konum;
      return Vec3d{static_cast<f64>(p[0]), static_cast<f64>(p[1]),
                   static_cast<f64>(p[2])};
    };
    c.ekle(v(m.indeksler[i]), v(m.indeksler[i + 1]), v(m.indeksler[i + 2]));
  }
  return c;
}

// Möller–Trumbore: keyfi `yon` ışını üçgeni ileride kesiyor mu (bariyer-sıkı).
static bool isinKesisim(Vec3d o, Vec3d yon, Vec3d a, Vec3d b, Vec3d cc) {
  const Vec3d e1 = b - a, e2 = cc - a;
  const Vec3d h = capraz(yon, e2);
  const f64 det = e1.nokta(h);
  if (std::abs(det) < 1e-12) return false;  // paralel
  const f64 inv = 1.0 / det;
  const Vec3d s = o - a;
  const f64 u = s.nokta(h) * inv;
  if (u < 0.0 || u > 1.0) return false;
  const Vec3d q = capraz(s, e1);
  const f64 vv = yon.nokta(q) * inv;
  if (vv < 0.0 || u + vv > 1.0) return false;
  const f64 t = e2.nokta(q) * inv;
  return t > 1e-9;  // yalnız ışın önü
}

bool noktaIcinde(const UcgenCorba& corba, Vec3d p) {
  return noktaIcindeHam(corba.vertler, p);
}

bool noktaIcindeHam(const std::vector<Vec3d>& V, Vec3d p) {
  if (V.size() < 3) return false;
  // Eksene-hizalı ışın kenar/köşe tekilliği → eğik yön (kutu yüzlerine paralel
  // değil). Sayaç kararsızsa (kenara isabet olasılığı) 3 farklı yön dener.
  static const Vec3d yonler[3] = {Vec3d{0.7139, 0.4113, 0.5667}.birim(),
                                  Vec3d{0.5231, 0.8112, 0.2593}.birim(),
                                  Vec3d{0.2887, 0.3491, 0.8913}.birim()};
  int icOy = 0;
  for (const Vec3d& yon : yonler) {
    u32 sayim = 0;
    for (size_t i = 0; i + 2 < V.size(); i += 3)
      if (isinKesisim(p, yon, V[i], V[i + 1], V[i + 2])) ++sayim;
    if (sayim & 1u) ++icOy;
  }
  return icOy >= 2;  // oy çokluğu (kenar tekilliğine dayanıklı)
}

f64 corbaHacim(const std::vector<Vec3d>& vertler) {
  // Diverjans teoremi: kapalı yüzey için 6V = Σ a·(b×c). Sarım yönü işareti
  // belirler → mutlak değer alınır (çorbanın yönü garanti değil).
  const size_t n = vertler.size() / 3;
  f64 alti = 0;
  for (size_t t = 0; t < n; ++t) {
    const Vec3d& a = vertler[3 * t];
    const Vec3d& b = vertler[3 * t + 1];
    const Vec3d& c = vertler[3 * t + 2];
    alti += a.nokta(b.capraz(c));
  }
  return std::abs(alti) / 6.0;
}

f64 kafesParcaKapsama(const VoxelKafes& kafes) {
  if (!kafes.gecerliMi() || kafes.parcaUcgen.empty()) return -1.0;
  const f64 gercek = corbaHacim(kafes.parcaUcgen);
  const f64 hm = kafes.hMesh();
  if (!(gercek > 0) || !(hm > 0) || !std::isfinite(gercek)) return -1.0;
  u32 ic = 0;
  for (auto t : kafes.tip)
    if (t == HucreTip::Ic) ++ic;
  return (static_cast<f64>(ic) * hm * hm * hm) / gercek;
}

f64 kafesMetalHacimM3(const VoxelKafes& kafes) {
  if (!kafes.gecerliMi()) return 0.0;
  u32 metal = 0;
  for (auto t : kafes.tip)
    if (t == HucreTip::Ic || t == HucreTip::Sinir) ++metal;
  return static_cast<f64>(metal) * kafes.hucreHacim();
}

bool& voxelSiniflaParalel() {
  static bool a = true;
  return a;
}

u32& voxelParalelEsik() {
  static u32 a = 4096;  // altında iş parçacığı kurulumu kazançtan pahalı
  return a;
}

void kafesiSinifla(VoxelKafes& kafes, const UcgenCorba& corba) {
  if (!kafes.gecerliMi()) return;
  const u32 nx = kafes.nx(), ny = kafes.ny(), nz = kafes.nz();
  // 1) İç/dış (ray parity, her hücre merkezi). Hücreler BİRBİRİNDEN BAĞIMSIZ →
  // dilim dilim paralel; her iş parçacığı yalnız kendi k aralığına yazar.
  const u32 isSayi =
      (voxelSiniflaParalel() && kafes.hucreSayi() >= voxelParalelEsik())
          ? std::max(1u, std::min(nz, std::thread::hardware_concurrency()))
          : 1u;
  const auto dilim = [&](u32 kBas, u32 kSon) {
    for (u32 k = kBas; k < kSon; ++k)
      for (u32 j = 0; j < ny; ++j)
        for (u32 i = 0; i < nx; ++i) {
          const u32 idx = kafes.indeks(i, j, k);
          kafes.tip[idx] =
              noktaIcinde(corba, kafes.merkez(i, j, k)) ? HucreTip::Ic : HucreTip::Dis;
        }
  };
  if (isSayi <= 1) {
    dilim(0, nz);
  } else {
    std::vector<std::thread> isler;
    isler.reserve(isSayi);
    const u32 pay = (nz + isSayi - 1) / isSayi;
    for (u32 t = 0; t < isSayi; ++t) {
      const u32 kBas = t * pay, kSon = std::min(nz, kBas + pay);
      if (kBas >= kSon) break;
      isler.emplace_back(dilim, kBas, kSon);
    }
    for (std::thread& i : isler) i.join();
  }
  // 2) Metale komşu dış hücre → Sinir (kalıp temas kabuğu).
  auto ic = [&](int i, int j, int k) -> bool {
    if (i < 0 || j < 0 || k < 0 || i >= (int)nx || j >= (int)ny || k >= (int)nz)
      return false;
    return kafes.tip[kafes.indeks(i, j, k)] == HucreTip::Ic;
  };
  for (u32 k = 0; k < nz; ++k)
    for (u32 j = 0; j < ny; ++j)
      for (u32 i = 0; i < nx; ++i) {
        const u32 idx = kafes.indeks(i, j, k);
        if (kafes.tip[idx] != HucreTip::Dis) continue;
        if (ic(i - 1, j, k) || ic(i + 1, j, k) || ic(i, j - 1, k) ||
            ic(i, j + 1, k) || ic(i, j, k - 1) || ic(i, j, k + 1))
          kafes.tip[idx] = HucreTip::Sinir;
      }
}

// ---- Nokta → üçgen EN YAKIN NOKTA (Ericson RTCD §5.1.5) -------------------
// Noktanın kendisini döndürür (yalnız mesafeyi değil): keskin köşede en-yakın
// nokta parçanın KÖŞE VERTEX'idir → oraya itilen tüm küp köşeleri aynı yere
// oturur → SİVRİ köşe. (SDF gradyanıyla itmek bunu YUVARLAR: trilineer aradeğer
// köşedeki koni kırığını yumuşatır → ölçüldü, keskinlik 1.0 → 0.9640.)
Vec3d noktaUcgenEnYakin(const Vec3d& p, const Vec3d& a, const Vec3d& b, const Vec3d& c) {
  const Vec3d ab = b - a, ac = c - a, ap = p - a;
  const f64 d1 = ab.nokta(ap), d2 = ac.nokta(ap);
  if (d1 <= 0 && d2 <= 0) return a;                           // köşe A
  const Vec3d bp = p - b;
  const f64 d3 = ab.nokta(bp), d4 = ac.nokta(bp);
  if (d3 >= 0 && d4 <= d3) return b;                          // köşe B
  const f64 vc = d1 * d4 - d3 * d2;
  if (vc <= 0 && d1 >= 0 && d3 <= 0) return a + ab * (d1 / (d1 - d3));  // kenar AB
  const Vec3d cp = p - c;
  const f64 d5 = ab.nokta(cp), d6 = ac.nokta(cp);
  if (d6 >= 0 && d5 <= d6) return c;                          // köşe C
  const f64 vb = d5 * d2 - d1 * d6;
  if (vb <= 0 && d2 >= 0 && d6 <= 0) return a + ac * (d2 / (d2 - d6));  // kenar AC
  const f64 va = d3 * d6 - d5 * d4;
  if (va <= 0 && (d4 - d3) >= 0 && (d5 - d6) >= 0)            // kenar BC
    return b + (c - b) * ((d4 - d3) / ((d4 - d3) + (d5 - d6)));
  const f64 denom = 1.0 / (va + vb + vc);                     // iç yüz
  return a + ab * (vb * denom) + ac * (vc * denom);
}

f64 noktaUcgenMesafe(const Vec3d& p, const Vec3d& a, const Vec3d& b, const Vec3d& c) {
  return (p - noktaUcgenEnYakin(p, a, b, c)).uzunluk();
}

Vec3d corbayaEnYakin(const UcgenCorba& corba, const Vec3d& p, f64 aramaYaricap) {
  const u32 n = corba.ucgenSayi();
  Vec3d enIyi = p;
  f64 enIyiD2 = 1e300;
  for (u32 t = 0; t < n; ++t) {
    const Vec3d& a = corba.vertler[3 * t];
    const Vec3d& b = corba.vertler[3 * t + 1];
    const Vec3d& c = corba.vertler[3 * t + 2];
    if (aramaYaricap > 0) {  // ucuz AABB ön-elemesi
      const f64 mnx = std::min({a.x, b.x, c.x}), mxx = std::max({a.x, b.x, c.x});
      if (p.x < mnx - aramaYaricap || p.x > mxx + aramaYaricap) continue;
      const f64 mny = std::min({a.y, b.y, c.y}), mxy = std::max({a.y, b.y, c.y});
      if (p.y < mny - aramaYaricap || p.y > mxy + aramaYaricap) continue;
      const f64 mnz = std::min({a.z, b.z, c.z}), mxz = std::max({a.z, b.z, c.z});
      if (p.z < mnz - aramaYaricap || p.z > mxz + aramaYaricap) continue;
    }
    const Vec3d q = noktaUcgenEnYakin(p, a, b, c);
    const f64 d2 = (p - q).uzunlukKare();
    if (d2 < enIyiD2) { enIyiD2 = d2; enIyi = q; }
  }
  return enIyiD2 >= 1e299 ? p : enIyi;  // hiç aday yok → dokunma
}

// ---- Gerçek parça SDF (yüzey kabuğunda exact, uzakta işaret) ---------------
void parcaSdfHesapla(VoxelKafes& kafes, const UcgenCorba& corba, f32 bantDisiDeger) {
  const u32 nx = kafes.nx(), ny = kafes.ny(), nz = kafes.nz();
  const u32 N = kafes.hucreSayi();
  const u32 nUcgen = corba.ucgenSayi();
  if (!kafes.gecerliMi() || kafes.tip.size() != N || nUcgen == 0) {
    kafes.parcaSdf.clear();
    return;
  }
  // Metal-lik maskesi = YALNIZ Ic. Sinir, parça yüzeyinin DIŞINDAKİ kalıp-temas
  // kabuğudur (bkz. sivi_yuzey parcaSdfKur notu) → maskeye alınırsa sıfır seviyesi
  // bir hücre dışarı kayar.
  auto metalMi = [&](u32 idx) { return kafes.tip[idx] == HucreTip::Ic; };

  // Kabuk: 6-komşusundan biri farklı metal-lik'te (grid kenarı da sınır sayılır).
  // MC yalnız bu hücreler arasındaki kenarlarda sıfır geçişi arar → exact mesafe
  // yalnız burada gerekir. Kabuk ≈ yüzey_alanı/h² (dev parçada bile küçük küme).
  std::vector<u8> kabuk(N, 0);
  for (u32 k = 0; k < nz; ++k)
    for (u32 j = 0; j < ny; ++j)
      for (u32 i = 0; i < nx; ++i) {
        const u32 idx = kafes.indeks(i, j, k);
        const bool m = metalMi(idx);
        auto farkli = [&](i64 a, i64 b, i64 c) {
          if (a < 0 || b < 0 || c < 0 || a >= nx || b >= ny || c >= nz) return m;
          return metalMi(kafes.indeks((u32)a, (u32)b, (u32)c)) != m;
        };
        if (farkli(i - 1, j, k) || farkli(i + 1, j, k) || farkli(i, j - 1, k) ||
            farkli(i, j + 1, k) || farkli(i, j, k - 1) || farkli(i, j, k + 1))
          kabuk[idx] = 1;
      }

  // Hücre boyutu MESH biriminde (merkez() mesh birimi döner; h() ise metre).
  const Vec3d m0 = kafes.merkez(0, 0, 0);
  const f64 hMesh = nx > 1 ? std::fabs(kafes.merkez(1, 0, 0).x - m0.x)
                           : std::cbrt(std::max(kafes.hucreHacim(), 1e-30));
  if (!(hMesh > 0)) { kafes.parcaSdf.clear(); return; }

  // Çorbayı sakla: görselleştirme küp köşelerini GERÇEK en-yakın-noktaya iter
  // (gradyan yerine) → keskin köşe korunur. Kopya maliyeti ~72 bayt/üçgen.
  kafes.parcaUcgen = corba.vertler;
  kafes.parcaSdf.assign(N, 0.0f);
  const f64 arama = 2.0 * hMesh;  // üçgen AABB ön-elemesi (kabuk hücresi çapı payı)
#pragma omp parallel for schedule(dynamic, 64) if (N > 4096)
  for (long long ii = 0; ii < (long long)N; ++ii) {
    const u32 idx = (u32)ii;
    const bool ic = metalMi(idx);
    if (!kabuk[idx]) {  // uzak alan: yalnız işaret önemli (MC sıfır geçişi aramaz)
      kafes.parcaSdf[idx] = ic ? -bantDisiDeger : bantDisiDeger;
      continue;
    }
    const u32 i = idx % nx, j = (idx / nx) % ny, k = idx / (nx * ny);
    const Vec3d p = kafes.merkez(i, j, k);
    f64 enYakin = 1e300;
    for (u32 t = 0; t < nUcgen; ++t) {
      const Vec3d& a = corba.vertler[3 * t];
      const Vec3d& b = corba.vertler[3 * t + 1];
      const Vec3d& c = corba.vertler[3 * t + 2];
      // Ucuz AABB ön-elemesi: üçgen bu hücreden arama yarıçapından uzaksa atla.
      const f64 mnx = std::min({a.x, b.x, c.x}), mxx = std::max({a.x, b.x, c.x});
      if (p.x < mnx - arama || p.x > mxx + arama) continue;
      const f64 mny = std::min({a.y, b.y, c.y}), mxy = std::max({a.y, b.y, c.y});
      if (p.y < mny - arama || p.y > mxy + arama) continue;
      const f64 mnz = std::min({a.z, b.z, c.z}), mxz = std::max({a.z, b.z, c.z});
      if (p.z < mnz - arama || p.z > mxz + arama) continue;
      enYakin = std::min(enYakin, noktaUcgenMesafe(p, a, b, c));
      if (enYakin <= 1e-12) break;
    }
    // Ön-eleme hepsini elediyse (kabuk hücresi ama yüzey uzak — ince boşluk vb.)
    // güvenli tarafa düş: bant dışı değer, işaret tip'ten.
    const f32 d = (enYakin >= 1e299) ? bantDisiDeger
                                     : (f32)(enYakin / hMesh);  // → hücre birimi
    kafes.parcaSdf[idx] = ic ? -d : d;
  }
}

}  // namespace fiko::dokum