#include <BRepFill_CompatibleWires.hxx>
#include <BRepBuilderAPI_MakeEdge.hxx>
#include <BRepBuilderAPI_MakeWire.hxx>
#include <BRepTools_WireExplorer.hxx>
#include <BRepAdaptor_Curve.hxx>
#include <BRepCheck_Analyzer.hxx>
#include <BRepTools.hxx>
#include <Geom_BezierCurve.hxx>
#include <TColgp_Array1OfPnt.hxx>
#include <TopTools_SequenceOfShape.hxx>
#include <TopoDS.hxx>
#include <TopoDS_Wire.hxx>
#include <Standard_Version.hxx>
#include <Standard_Failure.hxx>
#include <gp_Vec.hxx>
#include <algorithm>
#include <filesystem>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <string>
#include <vector>

using Samples = std::vector<gp_Pnt>;

const char* OrientationName(TopAbs_Orientation value)
{
  switch (value) {
    case TopAbs_FORWARD: return "FORWARD";
    case TopAbs_REVERSED: return "REVERSED";
    case TopAbs_INTERNAL: return "INTERNAL";
    case TopAbs_EXTERNAL: return "EXTERNAL";
  }
  return "UNKNOWN";
}

TopoDS_Wire WrapEdge(const TopoDS_Edge& edge)
{
  BRepBuilderAPI_MakeWire builder(edge);
  if (!builder.IsDone()) throw std::runtime_error("Cannot build wire");
  return builder.Wire();
}

int EdgeCount(const TopoDS_Wire& wire)
{
  int result = 0;
  for (BRepTools_WireExplorer it(wire); it.More(); it.Next()) ++result;
  return result;
}

// This demo samples by normalized parameter for drawing, not by arc length.
Samples SampleTraversal(const TopoDS_Wire& wire)
{
  Samples result;
  for (BRepTools_WireExplorer it(wire); it.More(); it.Next()) {
    const TopoDS_Edge edge = it.Current();
    BRepAdaptor_Curve curve(edge); // Includes the edge location.
    for (int i = 0; i <= 200; ++i) {
      double fraction = double(i) / 200.0;
      if (edge.Orientation() == TopAbs_REVERSED) fraction = 1.0 - fraction;
      const double u = curve.FirstParameter()
        + fraction * (curve.LastParameter() - curve.FirstParameter());
      result.push_back(curve.Value(u));
    }
  }
  return result;
}

void WritePoints(std::ostream& os, const Samples& points)
{
  os << '[';
  for (std::size_t i = 0; i < points.size(); ++i) {
    if (i != 0) os << ',';
    os << '[' << points[i].X() << ',' << points[i].Y() << ',' << points[i].Z() << ']';
  }
  os << ']';
}

// Finite sampled-set comparison, not an exact curve Hausdorff distance.
double SampledSetDifference(const Samples& a, const Samples& b)
{
  auto directed = [](const Samples& source, const Samples& target) {
    double maximum = 0.0;
    for (const gp_Pnt& p : source) {
      double minimum = 1e100;
      for (const gp_Pnt& q : target) minimum = std::min(minimum, p.Distance(q));
      maximum = std::max(maximum, minimum);
    }
    return maximum;
  };
  return std::max(directed(a, b), directed(b, a));
}

int Run(const std::filesystem::path& directory)
{
  std::filesystem::create_directories(directory);
  const auto path = [&](const char* name) { return (directory / name).string(); };

  // Two translated, non-planar cubic Bezier curves with identical parameterization.
  TColgp_Array1OfPnt polesA(1, 4), polesB(1, 4);
  polesA(1) = gp_Pnt(-5.0, 0.0, 0.0);
  polesA(2) = gp_Pnt(-2.5, 4.2, 1.8);
  polesA(3) = gp_Pnt( 2.5,-4.2,-1.2);
  polesA(4) = gp_Pnt( 5.0, 0.0, 0.8);
  for (int i = 1; i <= 4; ++i)
    polesB(i) = polesA(i).Translated(gp_Vec(0.0, 4.6, 3.0));

  Handle(Geom_BezierCurve) curveA = new Geom_BezierCurve(polesA);
  Handle(Geom_BezierCurve) curveB = new Geom_BezierCurve(polesB);
  BRepBuilderAPI_MakeEdge makeA(curveA), makeB(curveB);
  if (!makeA.IsDone() || !makeB.IsDone()) throw std::runtime_error("Cannot build edges");
  TopoDS_Edge edgeA = makeA.Edge(), edgeB = makeB.Edge();
  edgeB.Reverse(); // Prepare opposite traversal before calling the algorithm.
  const TopoDS_Wire wireA = WrapEdge(edgeA), wireB = WrapEdge(edgeB);

  TopTools_SequenceOfShape sections;
  sections.Append(wireA);
  sections.Append(wireB);
  BRepFill_CompatibleWires compatible(sections);
  compatible.Perform(); // No manual reversal is applied to the outputs.
  if (!compatible.IsDone()) {
    std::cerr << "Compatibility failed, status=" << int(compatible.GetStatus()) << '\n';
    return 2;
  }

  // Copy the returned sequence before the algorithm object is destroyed or reinitialized.
  const TopTools_SequenceOfShape adjusted = compatible.Shape();
  if (adjusted.Length() != 2) throw std::runtime_error("Unexpected number of wires");
  const TopoDS_Wire afterA = TopoDS::Wire(adjusted.Value(1));
  const TopoDS_Wire afterB = TopoDS::Wire(adjusted.Value(2));
  const Samples pa = SampleTraversal(wireA), pb = SampleTraversal(wireB);
  const Samples qa = SampleTraversal(afterA), qb = SampleTraversal(afterB);
  BRepTools_WireExplorer a(wireA), b(wireB), aa(afterA), bb(afterB);
  const bool sameEdgeTShapes = a.Current().IsSame(aa.Current())
    && b.Current().IsSame(bb.Current());
  const bool valid = BRepCheck_Analyzer(afterA).IsValid()
    && BRepCheck_Analyzer(afterB).IsValid();
  const double beforeDot = gp_Vec(pa.front(), pa.back()).Normalized()
    .Dot(gp_Vec(pb.front(), pb.back()).Normalized());
  const double afterDot = gp_Vec(qa.front(), qa.back()).Normalized()
    .Dot(gp_Vec(qb.front(), qb.back()).Normalized());
  const double displacement = std::max(SampledSetDifference(pa, qa), SampledSetDifference(pb, qb));

  std::ofstream json(path("curve-orientation-data.json"));
  if (!json) throw std::runtime_error("Cannot open JSON output");
  json << std::setprecision(12) << std::boolalpha;
  json << "{\"occtVersion\":\"" << OCC_VERSION_COMPLETE
       << "\",\"algorithm\":\"BRepFill_CompatibleWires::Perform\",\"isDone\":true"
       << ",\"validOutput\":" << valid << ",\"sameEdgeTShapes\":" << sameEdgeTShapes
       << ",\"beforeEndToEndDot\":" << beforeDot << ",\"afterEndToEndDot\":" << afterDot
       << ",\"sampledGeometryDisplacement\":" << displacement
       << ",\"inputEdgeCounts\":[" << EdgeCount(wireA) << ',' << EdgeCount(wireB)
       << "],\"outputEdgeCounts\":[" << EdgeCount(afterA) << ',' << EdgeCount(afterB)
       << "],\"before\":{\"a\":";
  WritePoints(json, pa); json << ",\"b\":"; WritePoints(json, pb);
  json << "},\"after\":{\"a\":";
  WritePoints(json, qa); json << ",\"b\":"; WritePoints(json, qb); json << "}}\n";

  BRepTools::Write(wireA, path("edge-a-before.brep").c_str());
  BRepTools::Write(wireB, path("edge-b-before.brep").c_str());
  BRepTools::Write(afterA, path("edge-a-after.brep").c_str());
  BRepTools::Write(afterB, path("edge-b-after.brep").c_str());

  std::cout << std::boolalpha << "OCCT " << OCC_VERSION_COMPLETE
    << "\nIsDone=" << compatible.IsDone() << ", valid=" << valid
    << "\nA " << OrientationName(a.Current().Orientation()) << " -> " << OrientationName(aa.Current().Orientation())
    << "\nB " << OrientationName(b.Current().Orientation()) << " -> " << OrientationName(bb.Current().Orientation())
    << "\nSameEdgeTShapes=" << sameEdgeTShapes
    << "\nEdges A " << EdgeCount(wireA) << " -> " << EdgeCount(afterA)
    << "; B " << EdgeCount(wireB) << " -> " << EdgeCount(afterB)
    << "\nDemo chord dot=" << beforeDot << " -> " << afterDot
    << "\nSampled geometry displacement=" << displacement << '\n';
  if (!sameEdgeTShapes || !valid || beforeDot > -0.99 || afterDot < 0.99 || displacement > 1e-9)
    return 3;
  return 0;
}

int main(int argc, char** argv)
{
  try {
    return Run(argc > 1 ? std::filesystem::path(argv[1]) : std::filesystem::path("run-results"));
  } catch (const Standard_Failure& failure) {
    std::cerr << "OCCT exception: " << failure.GetMessageString() << '\n';
  } catch (const std::exception& failure) {
    std::cerr << "Exception: " << failure.what() << '\n';
  }
  return 1;
}