Seminar5

// Geometry3D.hpp

#pragma once

#include <cmath>
#include <iostream>

using ld = long double;

using namespace std;

struct Point { ld x; ld y; ld z; };
struct Vector {
  ld x; ld y; ld z;
  ld Len() const { return sqrt(x * x + y * y + z * z); }
  ld Len2() const { return x * x + y * y + z * z; }
  void Normalize() { ld len = Len(); x /= len; y /= len; z /= len; }
};

Point operator+(const Point& p, const Vector& v) { return {p.x + v.x, p.y + v.y, p.z + v.z}; }
Point operator+(const Vector& v, const Point& p) { return {p.x + v.x, p.y + v.y, p.z + v.z}; }
Point operator-(const Point& p, const Vector& v) { return {p.x - v.x, p.y - v.y, p.z - v.z}; }
Vector operator-(const Point& p1, const Point& p2) { return {p1.x - p2.x, p1.y - p2.y, p1.z - p2.z}; }
Vector operator+(const Vector& v1, const Vector& v2) { return {v1.x + v2.x, v1.y + v2.y, v1.z + v2.z};}
Vector operator-(const Vector& v1, const Vector& v2) { return {v1.x - v2.x, v1.y - v2.y, v1.z - v2.z}; }
Vector operator*(const Vector& v, ld k) { return {k * v.x, k * v.y, k * v.z}; }
Vector operator*(ld k, const Vector& v) { return {v.x * k, v.y * k, v.z * k}; }
Vector operator/(const Vector& v, ld k) { return {v.x / k, v.y / k, v.z / k}; }
ld DotProduct(const Vector& v1, const Vector& v2) { return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z; }
Vector CrossProduct(const Vector& v1, const Vector& v2) {
  return {v1.y * v2.z - v1.z * v2.y, v1.z * v2.x - v1.x * v2.z, v1.x * v2.y - v1.y * v2.x};
}
istream& operator>>(istream& in, Point& v) { in >> v.x >> v.y >> v.z; return in; }
istream& operator>>(istream& in, Vector& v) { in >> v.x >> v.y >> v.z; return in; }
ostream& operator<<(ostream& out, const Point& v) { out << v.x << " " << v.y << " " << v.z; return out; }
ostream& operator<<(ostream& out, const Vector& v) { out << v.x << " " << v.y << " " << v.z; return out; }


// --------------------------------------------------------------------------------------------------------
// DCEL.hpp

#pragma once

#include "Geometry3D.hpp"

#include <vector>
#include <memory>
#include <functional>

using namespace std;

// Forward declarations
struct HalfEdge;
struct Vertex;
struct Face;

// Vertex structure - represents a point in the DCEL
struct Vertex : public Point {
  HalfEdge* incident_edge = nullptr;  // One of the half-edges originating from this vertex
  int index = -1;  // Index in the vertices array

  Vertex() = default;
  Vertex(ld x_, ld y_, ld z_) { x = x_; y = y_; z = z_; }
  Vertex(const Point& p) { x = p.x; y = p.y; z = p.z; }
};

// HalfEdge structure - the main component of DCEL
struct HalfEdge {
  Vertex* origin = nullptr;      // Origin vertex of this half-edge
  HalfEdge* twin = nullptr;      // Twin half-edge (opposite direction)
  HalfEdge* next = nullptr;      // Next half-edge in the face (counter-clockwise)
  HalfEdge* prev = nullptr;      // Previous half-edge in the face
  Face* incident_face = nullptr; // Face to the left of this half-edge
  int index = -1;                // Index in the edges array

  // Get the destination vertex
  Vertex* destination() const {
    return twin ? twin->origin : nullptr;
  }
};

// Face structure - represents a polygon face
struct Face {
  HalfEdge* outer_edge = nullptr;  // One of the half-edges on the outer boundary
  vector<HalfEdge*> inner_edges;   // Half-edges on inner boundaries (holes)
  int index = -1;                  // Index in the faces array
  bool is_outer = false;           // True if this is the unbounded outer face

  // Calculate the number of edges in this face
  int EdgeCount() const {
    if (!outer_edge) return 0;
    int count = 0;
    HalfEdge* e = outer_edge;
    do {
      count++;
      e = e->next;
    } while (e && e != outer_edge);
    return count;
  }

  // Get all vertices of this face
  vector<Vertex*> GetVertices() const {
    vector<Vertex*> result;
    if (!outer_edge) return result;
    HalfEdge* e = outer_edge;
    do {
      result.push_back(e->origin);
      e = e->next;
    } while (e && e != outer_edge);
    return result;
  }

  // Calculate the normal vector of the face (assuming planar face)
  Vector Normal() const {
    if (!outer_edge || !outer_edge->next) return {0, 0, 0};
    Vertex* v0 = outer_edge->origin;
    Vertex* v1 = outer_edge->destination();
    Vertex* v2 = outer_edge->next->destination();
    if (!v0 || !v1 || !v2) return {0, 0, 0};
    Vector e1 = *v1 - *v0;
    Vector e2 = *v2 - *v0;
    return CrossProduct(e1, e2);
  }
};

// DCEL class - Doubly Connected Edge List
class DCEL {
 public:
  DCEL() = default;

  // Accessors
  const vector<unique_ptr<HalfEdge>>& Edges() const { return edges_; }
  const vector<unique_ptr<Vertex>>& Vertices() const { return vertices_; }
  const vector<unique_ptr<Face>>& Faces() const { return faces_; }

  size_t NumEdges() const { return edges_.size(); }
  size_t NumVertices() const { return vertices_.size(); }
  size_t NumFaces() const { return faces_.size(); }

  // Create a new vertex
  Vertex* CreateVertex(ld x, ld y, ld z) {
    auto v = make_unique<Vertex>(x, y, z);
    v->index = static_cast<int>(vertices_.size());
    Vertex* ptr = v.get();
    vertices_.push_back(std::move(v));
    return ptr;
  }

  Vertex* CreateVertex(const Point& p) {
    return CreateVertex(p.x, p.y, p.z);
  }

  // Create a new half-edge pair (edge and its twin)
  pair<HalfEdge*, HalfEdge*> CreateEdgePair() {
    auto e1 = make_unique<HalfEdge>();
    auto e2 = make_unique<HalfEdge>();

    e1->index = static_cast<int>(edges_.size());
    e2->index = static_cast<int>(edges_.size() + 1);

    e1->twin = e2.get();
    e2->twin = e1.get();

    HalfEdge* ptr1 = e1.get();
    HalfEdge* ptr2 = e2.get();

    edges_.push_back(std::move(e1));
    edges_.push_back(std::move(e2));

    return {ptr1, ptr2};
  }

  // Create a new face
  Face* CreateFace() {
    auto f = make_unique<Face>();
    f->index = static_cast<int>(faces_.size());
    Face* ptr = f.get();
    faces_.push_back(std::move(f));
    return ptr;
  }

  // Connect half-edges in sequence (sets next/prev pointers)
  static void ConnectEdges(HalfEdge* e1, HalfEdge* e2) {
    if (e1) e1->next = e2;
    if (e2) e2->prev = e1;
  }

  // Set the origin of a half-edge and update vertex's incident edge if needed
  static void SetOrigin(HalfEdge* edge, Vertex* vertex) {
    if (!edge) return;
    edge->origin = vertex;
    if (vertex && !vertex->incident_edge) {
      vertex->incident_edge = edge;
    }
  }

  // Set the face for a half-edge and update face's outer edge if needed
  static void SetFace(HalfEdge* edge, Face* face) {
    if (!edge) return;
    edge->incident_face = face;
    if (face && !face->outer_edge) {
      face->outer_edge = edge;
    }
  }

  // Set face for all edges in a cycle starting from edge
  static void SetFaceForCycle(HalfEdge* start, Face* face) {
    if (!start || !face) return;
    HalfEdge* e = start;
    do {
      e->incident_face = face;
      e = e->next;
    } while (e && e != start);
    face->outer_edge = start;
  }

  // Build a triangular face from three vertices
  Face* BuildTriangle(Vertex* v0, Vertex* v1, Vertex* v2) {
    auto [e01, e10] = CreateEdgePair();
    auto [e12, e21] = CreateEdgePair();
    auto [e20, e02] = CreateEdgePair();

    // Set origins
    SetOrigin(e01, v0);
    SetOrigin(e10, v1);
    SetOrigin(e12, v1);
    SetOrigin(e21, v2);
    SetOrigin(e20, v2);
    SetOrigin(e02, v0);

    // Connect edges for the inner face (counter-clockwise: v0 -> v1 -> v2)
    ConnectEdges(e01, e12);
    ConnectEdges(e12, e20);
    ConnectEdges(e20, e01);

    // Create and set the face
    Face* face = CreateFace();
    SetFaceForCycle(e01, face);

    return face;
  }

  // Build a polygon face from a list of vertices
  Face* BuildPolygon(const vector<Vertex*>& vertices) {
    if (vertices.size() < 3) return nullptr;

    size_t n = vertices.size();
    vector<HalfEdge*> inner_edges(n);
    vector<HalfEdge*> outer_edges(n);

    // Create all edge pairs
    for (size_t i = 0; i < n; ++i) {
      auto [inner, outer] = CreateEdgePair();
      inner_edges[i] = inner;
      outer_edges[i] = outer;
    }

    // Set origins and connect edges
    for (size_t i = 0; i < n; ++i) {
      size_t next_i = (i + 1) % n;

      SetOrigin(inner_edges[i], vertices[i]);
      SetOrigin(outer_edges[i], vertices[next_i]);

      ConnectEdges(inner_edges[i], inner_edges[next_i]);
      ConnectEdges(outer_edges[next_i], outer_edges[i]);
    }

    // Create and set the face
    Face* face = CreateFace();
    SetFaceForCycle(inner_edges[0], face);

    return face;
  }

  // Iterate over all half-edges around a vertex (outgoing edges)
  void ForEachEdgeAroundVertex(Vertex* v, const function<void(HalfEdge*)>& func) const {
    if (!v || !v->incident_edge) return;
    HalfEdge* start = v->incident_edge;
    HalfEdge* e = start;
    do {
      func(e);
      if (!e->twin) break;
      e = e->twin->next;
    } while (e && e != start);
  }

  // Iterate over all half-edges of a face
  void ForEachEdgeOfFace(Face* f, const function<void(HalfEdge*)>& func) const {
    if (!f || !f->outer_edge) return;
    HalfEdge* start = f->outer_edge;
    HalfEdge* e = start;
    do {
      func(e);
      e = e->next;
    } while (e && e != start);
  }

  // Get all faces adjacent to a vertex
  vector<Face*> GetAdjacentFaces(Vertex* v) const {
    vector<Face*> result;
    ForEachEdgeAroundVertex(v, [&result](HalfEdge* e) {
      if (e->incident_face && !e->incident_face->is_outer) {
        result.push_back(e->incident_face);
      }
    });
    return result;
  }

  // Get all vertices adjacent to a vertex
  vector<Vertex*> GetAdjacentVertices(Vertex* v) const {
    vector<Vertex*> result;
    ForEachEdgeAroundVertex(v, [&result](HalfEdge* e) {
      if (e->destination()) {
        result.push_back(e->destination());
      }
    });
    return result;
  }

  // Check if the DCEL is valid (basic consistency checks)
  bool IsValid() const {
    for (const auto& e : edges_) {
      // Check twin relationship
      if (!e->twin || e->twin->twin != e.get()) {
        return false;
      }
      // Check next/prev relationship
      if (e->next && e->next->prev != e.get()) {
        return false;
      }
      if (e->prev && e->prev->next != e.get()) {
        return false;
      }
    }
    return true;
  }

  // Clear all data
  void Clear() {
    edges_.clear();
    vertices_.clear();
    faces_.clear();
  }

 private:
  vector<unique_ptr<HalfEdge>> edges_;
  vector<unique_ptr<Vertex>> vertices_;
  vector<unique_ptr<Face>> faces_;
};

// ----------------------------------------------------------------------------------------------------------------------
// main.cpp

#include <algorithm>
#include <array>
#include <cassert>
#include <climits>
#include <cmath>
#include <cstring>
#include <iomanip>
#include <iostream>
#include <limits>
#include <list>
#include <map>
#include <numeric>
#include <queue>
#include <random>
#include <set>
#include <sstream>
#include <stack>
#include <string>
#include <tuple>
#include <unordered_map>
#include <unordered_set>
#include <vector>

#include "DCEL.hpp"

#define all(x) x.begin(), x.end()
using namespace std;
using ll = long long;
using ld = long double;
const double eps = 1e-10;

// Test helper function
void TestPassed(const string& test_name) {
  cout << "[PASSED] " << test_name << endl;
}

void TestFailed(const string& test_name, const string& reason) {
  cout << "[FAILED] " << test_name << ": " << reason << endl;
}

// Test 1: Basic vertex creation
void TestVertexCreation() {
  DCEL dcel;

  Vertex* v1 = dcel.CreateVertex(0.0, 0.0, 0.0);
  Vertex* v2 = dcel.CreateVertex(1.0, 0.0, 0.0);
  Vertex* v3 = dcel.CreateVertex(0.5, 1.0, 0.0);

  assert(dcel.NumVertices() == 3);
  assert(v1->index == 0);
  assert(v2->index == 1);
  assert(v3->index == 2);
  assert(v1->x == 0.0 && v1->y == 0.0 && v1->z == 0.0);
  assert(v2->x == 1.0 && v2->y == 0.0 && v2->z == 0.0);
  assert(v3->x == 0.5 && v3->y == 1.0 && v3->z == 0.0);

  TestPassed("TestVertexCreation");
}

// Test 2: Edge pair creation
void TestEdgePairCreation() {
  DCEL dcel;

  auto [e1, e2] = dcel.CreateEdgePair();

  assert(dcel.NumEdges() == 2);
  assert(e1->twin == e2);
  assert(e2->twin == e1);
  assert(e1->index == 0);
  assert(e2->index == 1);

  TestPassed("TestEdgePairCreation");
}

// Test 3: Face creation
void TestFaceCreation() {
  DCEL dcel;

  Face* f = dcel.CreateFace();

  assert(dcel.NumFaces() == 1);
  assert(f->index == 0);
  assert(f->outer_edge == nullptr);

  TestPassed("TestFaceCreation");
}

// Test 4: Building a triangle
void TestBuildTriangle() {
  DCEL dcel;

  Vertex* v0 = dcel.CreateVertex(0.0, 0.0, 0.0);
  Vertex* v1 = dcel.CreateVertex(1.0, 0.0, 0.0);
  Vertex* v2 = dcel.CreateVertex(0.5, 1.0, 0.0);

  Face* face = dcel.BuildTriangle(v0, v1, v2);

  assert(face != nullptr);
  assert(face->outer_edge != nullptr);
  assert(face->EdgeCount() == 3);

  // Check that all vertices are in the face
  vector<Vertex*> face_vertices = face->GetVertices();
  assert(face_vertices.size() == 3);

  // Check that the cycle is closed
  HalfEdge* start = face->outer_edge;
  HalfEdge* e = start;
  int count = 0;
  do {
    assert(e->incident_face == face);
    assert(e->next != nullptr);
    assert(e->prev != nullptr);
    assert(e->next->prev == e);
    assert(e->prev->next == e);
    e = e->next;
    count++;
  } while (e != start && count < 10);
  assert(count == 3);

  TestPassed("TestBuildTriangle");
}

// Test 5: Building a polygon (square)
void TestBuildPolygon() {
  DCEL dcel;

  Vertex* v0 = dcel.CreateVertex(0.0, 0.0, 0.0);
  Vertex* v1 = dcel.CreateVertex(1.0, 0.0, 0.0);
  Vertex* v2 = dcel.CreateVertex(1.0, 1.0, 0.0);
  Vertex* v3 = dcel.CreateVertex(0.0, 1.0, 0.0);

  Face* face = dcel.BuildPolygon({v0, v1, v2, v3});

  assert(face != nullptr);
  assert(face->EdgeCount() == 4);

  vector<Vertex*> face_vertices = face->GetVertices();
  assert(face_vertices.size() == 4);

  TestPassed("TestBuildPolygon");
}

// Test 6: Edge connectivity
void TestEdgeConnectivity() {
  DCEL dcel;

  auto [e1, e2] = dcel.CreateEdgePair();
  auto [e3, e4] = dcel.CreateEdgePair();

  DCEL::ConnectEdges(e1, e3);

  assert(e1->next == e3);
  assert(e3->prev == e1);

  TestPassed("TestEdgeConnectivity");
}

// Test 7: Vertex incident edge
void TestVertexIncidentEdge() {
  DCEL dcel;

  Vertex* v = dcel.CreateVertex(0.0, 0.0, 0.0);
  auto [e1, e2] = dcel.CreateEdgePair();

  DCEL::SetOrigin(e1, v);

  assert(e1->origin == v);
  assert(v->incident_edge == e1);

  TestPassed("TestVertexIncidentEdge");
}

// Test 8: Face normal calculation
void TestFaceNormal() {
  DCEL dcel;

  // Create a triangle in the XY plane
  Vertex* v0 = dcel.CreateVertex(0.0, 0.0, 0.0);
  Vertex* v1 = dcel.CreateVertex(1.0, 0.0, 0.0);
  Vertex* v2 = dcel.CreateVertex(0.0, 1.0, 0.0);

  Face* face = dcel.BuildTriangle(v0, v1, v2);
  Vector normal = face->Normal();

  // Normal should point in Z direction (positive or negative depending on winding)
  assert(abs(normal.x) < eps);
  assert(abs(normal.y) < eps);
  assert(abs(normal.z) > eps);

  TestPassed("TestFaceNormal");
}

// Test 9: DCEL validity check
void TestDCELValidity() {
  DCEL dcel;

  Vertex* v0 = dcel.CreateVertex(0.0, 0.0, 0.0);
  Vertex* v1 = dcel.CreateVertex(1.0, 0.0, 0.0);
  Vertex* v2 = dcel.CreateVertex(0.5, 1.0, 0.0);

  dcel.BuildTriangle(v0, v1, v2);

  assert(dcel.IsValid());

  TestPassed("TestDCELValidity");
}

// Test 10: Destination vertex
void TestDestinationVertex() {
  DCEL dcel;

  Vertex* v0 = dcel.CreateVertex(0.0, 0.0, 0.0);
  Vertex* v1 = dcel.CreateVertex(1.0, 0.0, 0.0);

  auto [e1, e2] = dcel.CreateEdgePair();
  DCEL::SetOrigin(e1, v0);
  DCEL::SetOrigin(e2, v1);

  assert(e1->destination() == v1);
  assert(e2->destination() == v0);

  TestPassed("TestDestinationVertex");
}

// Test 11: Adjacent vertices
void TestAdjacentVertices() {
  DCEL dcel;

  Vertex* v0 = dcel.CreateVertex(0.0, 0.0, 0.0);
  Vertex* v1 = dcel.CreateVertex(1.0, 0.0, 0.0);
  Vertex* v2 = dcel.CreateVertex(0.5, 1.0, 0.0);

  dcel.BuildTriangle(v0, v1, v2);

  // v0 should be adjacent to v1 and v2
  vector<Vertex*> adjacent = dcel.GetAdjacentVertices(v0);
  assert(adjacent.size() >= 1);  // At least one adjacent vertex

  TestPassed("TestAdjacentVertices");
}

// Test 12: Clear DCEL
void TestClearDCEL() {
  DCEL dcel;

  dcel.CreateVertex(0.0, 0.0, 0.0);
  dcel.CreateVertex(1.0, 0.0, 0.0);
  dcel.CreateEdgePair();
  dcel.CreateFace();

  assert(dcel.NumVertices() == 2);
  assert(dcel.NumEdges() == 2);
  assert(dcel.NumFaces() == 1);

  dcel.Clear();

  assert(dcel.NumVertices() == 0);
  assert(dcel.NumEdges() == 0);
  assert(dcel.NumFaces() == 0);

  TestPassed("TestClearDCEL");
}

// Test 13: ForEachEdgeOfFace
void TestForEachEdgeOfFace() {
  DCEL dcel;

  Vertex* v0 = dcel.CreateVertex(0.0, 0.0, 0.0);
  Vertex* v1 = dcel.CreateVertex(1.0, 0.0, 0.0);
  Vertex* v2 = dcel.CreateVertex(0.5, 1.0, 0.0);

  Face* face = dcel.BuildTriangle(v0, v1, v2);

  int edge_count = 0;
  dcel.ForEachEdgeOfFace(face, [&edge_count](HalfEdge* e) {
    edge_count++;
    assert(e != nullptr);
  });

  assert(edge_count == 3);

  TestPassed("TestForEachEdgeOfFace");
}

// Test 14: Pentagon polygon
void TestPentagon() {
  DCEL dcel;

  // Create a regular pentagon
  vector<Vertex*> vertices;
  for (int i = 0; i < 5; ++i) {
    ld angle = 2.0 * M_PI * i / 5.0;
    vertices.push_back(dcel.CreateVertex(cos(angle), sin(angle), 0.0));
  }

  Face* face = dcel.BuildPolygon(vertices);

  assert(face != nullptr);
  assert(face->EdgeCount() == 5);
  assert(dcel.IsValid());

  TestPassed("TestPentagon");
}

// Test 15: Multiple triangles
void TestMultipleTriangles() {
  DCEL dcel;

  // Create two separate triangles
  Vertex* v0 = dcel.CreateVertex(0.0, 0.0, 0.0);
  Vertex* v1 = dcel.CreateVertex(1.0, 0.0, 0.0);
  Vertex* v2 = dcel.CreateVertex(0.5, 1.0, 0.0);

  Vertex* v3 = dcel.CreateVertex(2.0, 0.0, 0.0);
  Vertex* v4 = dcel.CreateVertex(3.0, 0.0, 0.0);
  Vertex* v5 = dcel.CreateVertex(2.5, 1.0, 0.0);

  Face* face1 = dcel.BuildTriangle(v0, v1, v2);
  Face* face2 = dcel.BuildTriangle(v3, v4, v5);

  assert(face1 != nullptr);
  assert(face2 != nullptr);
  assert(dcel.NumFaces() == 2);
  assert(dcel.NumVertices() == 6);
  assert(dcel.IsValid());

  TestPassed("TestMultipleTriangles");
}

void RunAllTests() {
  cout << "Running DCEL tests..." << endl;
  cout << "========================" << endl;

  TestVertexCreation();
  TestEdgePairCreation();
  TestFaceCreation();
  TestBuildTriangle();
  TestBuildPolygon();
  TestEdgeConnectivity();
  TestVertexIncidentEdge();
  TestFaceNormal();
  TestDCELValidity();
  TestDestinationVertex();
  TestAdjacentVertices();
  TestClearDCEL();
  TestForEachEdgeOfFace();
  TestPentagon();
  TestMultipleTriangles();

  cout << "========================" << endl;
  cout << "All tests passed!" << endl;
}

void Solve() {
  RunAllTests();
}

int main() {
  ios::sync_with_stdio(false);
  cin.tie(nullptr);
  cout.tie(nullptr);
  Solve();
  return 0;
}