#include <asiAlgo_AAG.h>
#include <asiAlgo_FeatureAttrAdjacency.h>
#include <asiAlgo_FeatureAttrAngle.h>
#include <asiAlgo_FeatureFaces.h>
#include <asiAlgo_RecognizeDrillHoles.h>
#include <asiAlgo_STEP.h>

#include <BRepAdaptor_Surface.hxx>
#include <BRepCheck_Analyzer.hxx>
#include <GeomAbs_SurfaceType.hxx>
#include <Precision.hxx>
#include <TColStd_MapIteratorOfPackedMapOfInteger.hxx>
#include <TopAbs_ShapeEnum.hxx>
#include <TopExp.hxx>
#include <TopTools_IndexedMapOfShape.hxx>
#include <gp_Ax1.hxx>
#include <gp_Cylinder.hxx>
#include <gp_Dir.hxx>
#include <gp_Pln.hxx>

#include <cmath>
#include <cstdlib>
#include <iostream>
#include <map>
#include <string>
#include <vector>

namespace
{
void printFeature(const asiAlgo_Feature& feature)
{
  std::cout << "[";
  bool first = true;
  for (TColStd_MapIteratorOfPackedMapOfInteger it(feature); it.More(); it.Next())
  {
    if (!first)
      std::cout << ", ";
    first = false;
    std::cout << it.Key();
  }
  std::cout << "]";
}

void printPackedMap(const TColStd_PackedMapOfInteger& values)
{
  std::cout << "[";
  bool first = true;
  for (TColStd_MapIteratorOfPackedMapOfInteger it(values); it.More(); it.Next())
  {
    if (!first)
      std::cout << ", ";
    first = false;
    std::cout << it.Key();
  }
  std::cout << "]";
}

const char* surfaceTypeName(const GeomAbs_SurfaceType type)
{
  switch (type)
  {
    case GeomAbs_Plane:
      return "plane";
    case GeomAbs_Cylinder:
      return "cylinder";
    case GeomAbs_Cone:
      return "cone";
    case GeomAbs_Sphere:
      return "sphere";
    case GeomAbs_Torus:
      return "torus";
    case GeomAbs_BezierSurface:
      return "bezier";
    case GeomAbs_BSplineSurface:
      return "bspline";
    case GeomAbs_SurfaceOfRevolution:
      return "revolution";
    case GeomAbs_SurfaceOfExtrusion:
      return "extrusion";
    case GeomAbs_OffsetSurface:
      return "offset";
    default:
      return "other";
  }
}

void printStringIntMap(const std::map<std::string, int>& values)
{
  std::cout << "{";
  bool first = true;
  for (const auto& item : values)
  {
    if (!first)
      std::cout << ", ";
    first = false;
    std::cout << "\"" << item.first << "\": " << item.second;
  }
  std::cout << "}";
}

void printGeometricRelation(
  const int faceId,
  const int neighborId,
  const std::string& type,
  const double residual,
  const char* source,
  const bool first)
{
  if (!first)
    std::cout << ",\n";
  std::cout << "    { \"faceIds\": [" << faceId << ", " << neighborId << "]"
            << ", \"type\": \"" << type << "\""
            << ", \"residual\": " << residual
            << ", \"source\": \"" << source << "\" }";
}
}

int main(int argc, char** argv)
{
  if (argc < 2)
  {
    std::cerr << "Usage: recognize_holes <step-file> [max-radius]\n";
    return 1;
  }

  const TCollection_AsciiString filename(argv[1]);
  const double radius = argc >= 3 ? std::atof(argv[2]) : Precision::Infinite();

  TopoDS_Shape shape;
  if (!asiAlgo_STEP::Import(filename, shape))
  {
    std::cerr << "Failed to import STEP: " << argv[1] << "\n";
    return 2;
  }

  TopTools_IndexedMapOfShape faces;
  TopExp::MapShapes(shape, TopAbs_FACE, faces);

  BRepCheck_Analyzer analyzer(shape);
  Handle(asiAlgo_AAG) aag = new asiAlgo_AAG(
    shape,
    false,
    1.0e-4,
    asiAlgo_AAG::CachedMap_All);

  asiAlgo_RecognizeDrillHoles recognizer(aag, true);
  if (!recognizer.Perform(radius))
  {
    std::cerr << "Hole recognition failed.\n";
    return 3;
  }

  const asiAlgo_Feature& holeFaceIds = recognizer.GetResultIndices();
  std::vector<asiAlgo_Feature> holes;
  aag->GetConnectedComponents(holeFaceIds, holes);

  std::map<std::string, int> surfaceSummary;
  std::map<std::string, int> angleSummary;
  std::map<std::string, int> geometricRelationSummary;
  for (int faceId = 1; faceId <= aag->GetNumberOfNodes(); ++faceId)
  {
    const TopoDS_Face& face = aag->GetFace(faceId);
    BRepAdaptor_Surface surface(face);
    surfaceSummary[surfaceTypeName(surface.GetType())]++;
  }

  std::cout << "{\n";
  std::cout << "  \"validBreP\": " << (analyzer.IsValid() ? "true" : "false") << ",\n";
  std::cout << "  \"faceCount\": " << faces.Extent() << ",\n";
  std::cout << "  \"aagNodeCount\": " << aag->GetNumberOfNodes() << ",\n";
  std::cout << "  \"holeFaceIds\": ";
  printFeature(holeFaceIds);
  std::cout << ",\n";
  std::cout << "  \"holeCount\": " << holes.size() << ",\n";
  std::cout << "  \"holes\": [\n";
  for (std::size_t i = 0; i < holes.size(); ++i)
  {
    std::cout << "    { \"index\": " << (i + 1) << ", \"faceIds\": ";
    printFeature(holes[i]);
    std::cout << " }";
    if (i + 1 < holes.size())
      std::cout << ",";
    std::cout << "\n";
  }
  std::cout << "  ],\n";

  std::cout << "  \"faces\": [\n";
  for (int faceId = 1; faceId <= aag->GetNumberOfNodes(); ++faceId)
  {
    const TopoDS_Face& face = aag->GetFace(faceId);
    BRepAdaptor_Surface surface(face);
    std::cout << "    { \"id\": " << faceId
              << ", \"surface\": \"" << surfaceTypeName(surface.GetType()) << "\""
              << ", \"neighbors\": ";
    if (aag->HasNeighbors(faceId))
      printFeature(aag->GetNeighbors(faceId));
    else
      std::cout << "[]";
    std::cout << " }";
    if (faceId < aag->GetNumberOfNodes())
      std::cout << ",";
    std::cout << "\n";
  }
  std::cout << "  ],\n";

  std::cout << "  \"adjacency\": [\n";
  bool firstAdjacency = true;
  std::vector<std::pair<int, int>> tangentPairs;
  for (int faceId = 1; faceId <= aag->GetNumberOfNodes(); ++faceId)
  {
    if (!aag->HasNeighbors(faceId))
      continue;
    const asiAlgo_Feature& neighbors = aag->GetNeighbors(faceId);
    for (TColStd_MapIteratorOfPackedMapOfInteger it(neighbors); it.More(); it.Next())
    {
      const int neighborId = it.Key();
      if (neighborId <= faceId)
        continue;

      const asiAlgo_AAG::t_arc arc(faceId, neighborId);
      Handle(asiAlgo_FeatureAttrAngle) angleAttr =
        Handle(asiAlgo_FeatureAttrAngle)::DownCast(aag->GetArcAttribute(arc));
      Handle(asiAlgo_FeatureAttrAdjacency) adjacencyAttr =
        Handle(asiAlgo_FeatureAttrAdjacency)::DownCast(aag->GetArcAttribute(arc));

      std::string angleType = "adjacent";
      double angleRad = 0.0;
      if (!angleAttr.IsNull())
      {
        angleType = asiAlgo_FeatureAngle::ToString(angleAttr->GetAngleType());
        angleRad = angleAttr->GetAngleRad();
      }
      angleSummary[angleType]++;
      if (angleType.find("smooth") != std::string::npos)
        tangentPairs.emplace_back(faceId, neighborId);

      if (!firstAdjacency)
        std::cout << ",\n";
      firstAdjacency = false;
      std::cout << "    { \"faceIds\": [" << faceId << ", " << neighborId << "]"
                << ", \"angleType\": \"" << angleType << "\""
                << ", \"angleRad\": " << angleRad
                << ", \"edgeIds\": ";
      if (!adjacencyAttr.IsNull())
        printPackedMap(adjacencyAttr->GetEdgeIndices());
      else
        std::cout << "[]";
      std::cout << " }";
    }
  }
  std::cout << "\n  ],\n";

  std::cout << "  \"geometricRelations\": [\n";
  bool firstRelation = true;
  const double angleTol = 1.0e-7;
  const double linearTol = 1.0e-4;
  const int geometricPairFaceLimit = 800;
  const bool runFullGeometricPairScan = aag->GetNumberOfNodes() <= geometricPairFaceLimit;
  if (runFullGeometricPairScan)
  {
    for (int faceId = 1; faceId <= aag->GetNumberOfNodes(); ++faceId)
    {
      const TopoDS_Face& leftFace = aag->GetFace(faceId);
      BRepAdaptor_Surface leftSurface(leftFace);
      for (int neighborId = faceId + 1; neighborId <= aag->GetNumberOfNodes(); ++neighborId)
      {
        const TopoDS_Face& rightFace = aag->GetFace(neighborId);
        BRepAdaptor_Surface rightSurface(rightFace);
        std::string relationType;
        double residual = 0.0;

        if (leftSurface.GetType() == GeomAbs_Plane && rightSurface.GetType() == GeomAbs_Plane)
        {
          const gp_Pln leftPlane = leftSurface.Plane();
          const gp_Pln rightPlane = rightSurface.Plane();
          const gp_Dir leftDir = leftPlane.Axis().Direction();
          const gp_Dir rightDir = rightPlane.Axis().Direction();
          if (leftDir.IsParallel(rightDir, angleTol))
          {
            residual = leftPlane.Distance(rightPlane.Location());
            relationType = residual <= linearTol ? "coplanar" : "parallel";
          }
          else if (leftDir.IsNormal(rightDir, angleTol))
          {
            residual = std::abs(leftDir.Dot(rightDir));
            relationType = "perpendicular";
          }
        }
        else if (leftSurface.GetType() == GeomAbs_Cylinder && rightSurface.GetType() == GeomAbs_Cylinder)
        {
          const gp_Cylinder leftCylinder = leftSurface.Cylinder();
          const gp_Cylinder rightCylinder = rightSurface.Cylinder();
          const gp_Ax1 leftAxis = leftCylinder.Axis();
          const gp_Ax1 rightAxis = rightCylinder.Axis();
          if (leftAxis.IsCoaxial(rightAxis, angleTol, linearTol) ||
              leftAxis.IsCoaxial(rightAxis.Reversed(), angleTol, linearTol))
          {
            residual = std::abs(leftCylinder.Radius() - rightCylinder.Radius());
            relationType = "coaxial";
          }
          else if (leftAxis.Direction().IsParallel(rightAxis.Direction(), angleTol))
          {
            residual = 0.0;
            relationType = "parallel_axis";
          }
        }

        if (!relationType.empty())
        {
          printGeometricRelation(faceId, neighborId, relationType, residual, "analysis-situs-probe", firstRelation);
          firstRelation = false;
          geometricRelationSummary[relationType]++;
        }
      }
    }
  }
  for (const auto& item : tangentPairs)
  {
    printGeometricRelation(item.first, item.second, "tangent", 0.0, "analysis-situs-aag-angle", firstRelation);
    firstRelation = false;
    geometricRelationSummary["tangent"]++;
  }
  std::cout << "\n  ],\n";

  std::cout << "  \"surfaceSummary\": ";
  printStringIntMap(surfaceSummary);
  std::cout << ",\n";
  std::cout << "  \"angleSummary\": ";
  printStringIntMap(angleSummary);
  std::cout << ",\n";
  std::cout << "  \"geometricRelationMode\": \""
            << (runFullGeometricPairScan ? "all-pairs" : "aag-smooth-only-large-model") << "\",\n";
  std::cout << "  \"geometricRelationSummary\": ";
  printStringIntMap(geometricRelationSummary);
  std::cout << "\n";
  std::cout << "}\n";

  return 0;
}