// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later

// primitives.hpp — 幾何プリミティブ (2D/3D)
//
// 移植元: lib++20 メディア/画像/図形/ (点, 直線, 線分, 平面, 球, 円柱, 円錐, 三角形, 四角形 等)
// 改良点:
//   - calx namespace, 英語 API
//   - std::array ベース (固定サイズ)
//   - constexpr 対応
//   - operator<=> (C++20)

#ifndef CALX_IMAGE_GEOMETRY_PRIMITIVES_HPP
#define CALX_IMAGE_GEOMETRY_PRIMITIVES_HPP

#include <array>
#include <cassert>
#include <cmath>
#include <ostream>

namespace calx::image {

    // ========================================================================
    // 2D プリミティブ
    // ========================================================================

    /// 2D 点
    template<typename T>
    struct Point2D {
        T x{}, y{};

        Point2D() = default;
        constexpr Point2D(T x, T y) : x(x), y(y) {}

        constexpr Point2D operator+(const Point2D& p) const { return {x + p.x, y + p.y}; }
        constexpr Point2D operator-(const Point2D& p) const { return {x - p.x, y - p.y}; }
        constexpr Point2D operator*(T s) const { return {x * s, y * s}; }
        friend constexpr Point2D operator*(T s, const Point2D& p) { return p * s; }

        bool operator==(const Point2D&) const = default;

        T norm() const { return std::sqrt(x * x + y * y); }
        T normSq() const { return x * x + y * y; }
    };

    /// 2 点間のユークリッド距離
    template<typename T>
    T distance(const Point2D<T>& a, const Point2D<T>& b) {
        return (b - a).norm();
    }

    /// 2D 内積
    template<typename T>
    constexpr T dot(const Point2D<T>& a, const Point2D<T>& b) {
        return a.x * b.x + a.y * b.y;
    }

    /// 2D 外積 (スカラー)
    template<typename T>
    constexpr T cross(const Point2D<T>& a, const Point2D<T>& b) {
        return a.x * b.y - a.y * b.x;
    }

    /// 2D 直線 (2 点で定義)
    template<typename T>
    struct Line2D {
        Point2D<T> p0, p1;

        Line2D() = default;
        constexpr Line2D(const Point2D<T>& a, const Point2D<T>& b) : p0(a), p1(b) {}

        /// 点 P から直線への距離
        T distanceTo(const Point2D<T>& p) const {
            auto v = p1 - p0;
            return std::abs(cross(p - p0, v)) / v.norm();
        }

        /// 直線上で P に最も近い点
        Point2D<T> closestPoint(const Point2D<T>& p) const {
            auto v = p1 - p0;
            T t = dot(p - p0, v) / dot(v, v);
            return p0 + t * v;
        }
    };

    /// 2D 線分
    template<typename T>
    struct Segment2D {
        Point2D<T> p0, p1;

        Segment2D() = default;
        constexpr Segment2D(const Point2D<T>& a, const Point2D<T>& b) : p0(a), p1(b) {}

        T length() const { return distance(p0, p1); }

        /// 点 P から線分への最短距離
        T distanceTo(const Point2D<T>& p) const {
            auto v = p1 - p0;
            auto w = p - p0;
            T c1 = dot(w, v);
            if (c1 <= T{0}) return distance(p, p0);
            T c2 = dot(v, v);
            if (c2 <= c1) return distance(p, p1);
            T t = c1 / c2;
            Point2D<T> proj = p0 + t * v;
            return distance(p, proj);
        }

        /// 2 線分の交差判定
        static bool intersects(const Segment2D& s1, const Segment2D& s2) {
            auto signedArea = [](const Point2D<T>& a, const Point2D<T>& b, const Point2D<T>& c) {
                return (a.x - c.x) * (b.y - c.y) - (b.x - c.x) * (a.y - c.y);
            };
            T d1 = signedArea(s1.p0, s1.p1, s2.p0);
            T d2 = signedArea(s1.p0, s1.p1, s2.p1);
            T d3 = signedArea(s2.p0, s2.p1, s1.p0);
            T d4 = signedArea(s2.p0, s2.p1, s1.p1);
            return d1 * d2 <= T{0} && d3 * d4 <= T{0};
        }
    };

    /// 2D 三角形
    template<typename T>
    struct Triangle2D {
        Point2D<T> p0, p1, p2;

        Triangle2D() = default;
        constexpr Triangle2D(const Point2D<T>& a, const Point2D<T>& b, const Point2D<T>& c)
            : p0(a), p1(b), p2(c) {}

        /// 面積 (符号付き)
        constexpr T signedArea() const {
            return ((p0.x - p2.x) * (p1.y - p2.y) - (p1.x - p2.x) * (p0.y - p2.y)) / T{2};
        }

        T area() const { return std::abs(signedArea()); }
    };

    // ========================================================================
    // 3D プリミティブ
    // ========================================================================

    /// 3D 点
    template<typename T>
    struct Point3D {
        T x{}, y{}, z{};

        Point3D() = default;
        constexpr Point3D(T x, T y, T z) : x(x), y(y), z(z) {}

        constexpr Point3D operator+(const Point3D& p) const { return {x + p.x, y + p.y, z + p.z}; }
        constexpr Point3D operator-(const Point3D& p) const { return {x - p.x, y - p.y, z - p.z}; }
        constexpr Point3D operator*(T s) const { return {x * s, y * s, z * s}; }
        friend constexpr Point3D operator*(T s, const Point3D& p) { return p * s; }

        bool operator==(const Point3D&) const = default;

        T norm() const { return std::sqrt(x * x + y * y + z * z); }
        T normSq() const { return x * x + y * y + z * z; }

        Point3D normalized() const {
            T n = norm();
            return {x / n, y / n, z / n};
        }
    };

    /// 3 点間のユークリッド距離
    template<typename T>
    T distance(const Point3D<T>& a, const Point3D<T>& b) {
        return (b - a).norm();
    }

    /// 3D 内積
    template<typename T>
    constexpr T dot(const Point3D<T>& a, const Point3D<T>& b) {
        return a.x * b.x + a.y * b.y + a.z * b.z;
    }

    /// 3D 外積
    template<typename T>
    constexpr Point3D<T> cross(const Point3D<T>& a, const Point3D<T>& 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};
    }

    /// 3D 直線 (2 点で定義)
    template<typename T>
    struct Line3D {
        Point3D<T> p0, p1;

        Line3D() = default;
        constexpr Line3D(const Point3D<T>& a, const Point3D<T>& b) : p0(a), p1(b) {}

        /// 点 P から直線への距離
        T distanceTo(const Point3D<T>& p) const {
            return cross(p - p0, p1 - p0).norm() / distance(p0, p1);
        }

        /// 直線上で P に最も近い点
        Point3D<T> closestPoint(const Point3D<T>& p) const {
            auto v = (p1 - p0).normalized();
            T d = dot(v, p - p0);
            return p0 + d * v;
        }
    };

    /// 3D 線分
    template<typename T>
    struct Segment3D {
        Point3D<T> p0, p1;

        Segment3D() = default;
        constexpr Segment3D(const Point3D<T>& a, const Point3D<T>& b) : p0(a), p1(b) {}

        T length() const { return distance(p0, p1); }
    };

    /// 3D 平面 (点 + 法線)
    template<typename T>
    struct Plane {
        Point3D<T> origin;
        Point3D<T> normal;

        Plane() = default;

        /// 法線と点から構成
        Plane(const Point3D<T>& o, const Point3D<T>& n) : origin(o), normal(n) {}

        /// 3 点から構成
        Plane(const Point3D<T>& a, const Point3D<T>& b, const Point3D<T>& c)
            : origin(a), normal(cross(b - a, c - a).normalized()) {}

        /// 点から平面への距離 (符号付き)
        T signedDistanceTo(const Point3D<T>& p) const {
            return dot(normal, p - origin);
        }

        T distanceTo(const Point3D<T>& p) const {
            return std::abs(signedDistanceTo(p));
        }

        /// 垂線の足
        Point3D<T> project(const Point3D<T>& p) const {
            T d = signedDistanceTo(p);
            return p - d * normal;
        }
    };

    /// 3D 三角形
    template<typename T>
    struct Triangle3D {
        Point3D<T> p0, p1, p2;

        Triangle3D() = default;
        constexpr Triangle3D(const Point3D<T>& a, const Point3D<T>& b, const Point3D<T>& c)
            : p0(a), p1(b), p2(c) {}

        /// 面積 (ヘロンの公式)
        T area() const {
            T a = distance(p0, p1);
            T b = distance(p1, p2);
            T c = distance(p2, p0);
            T s = (a + b + c) / T{2};
            return std::sqrt(s * (s - a) * (s - b) * (s - c));
        }

        /// 法線ベクトル (正規化済み)
        Point3D<T> normal() const {
            return cross(p1 - p0, p2 - p0).normalized();
        }
    };

    /// 3D 球
    template<typename T>
    struct Sphere {
        Point3D<T> center;
        T radius{};

        Sphere() = default;
        Sphere(const Point3D<T>& c, T r) : center(c), radius(r) {}

        bool contains(const Point3D<T>& p) const {
            return distance(center, p) <= radius;
        }
    };

    /// 3D 円柱 (2 端点 + 半径)
    template<typename T>
    struct Cylinder {
        Point3D<T> p0, p1;
        T radius{};

        Cylinder() = default;
        Cylinder(const Point3D<T>& a, const Point3D<T>& b, T r) : p0(a), p1(b), radius(r) {}

        T height() const { return distance(p0, p1); }
    };

    /// 3D 円錐 (底面中心, 頂点, 底面半径)
    template<typename T>
    struct Cone {
        Point3D<T> base;
        Point3D<T> apex;
        T radius{};

        Cone() = default;
        Cone(const Point3D<T>& b, const Point3D<T>& a, T r) : base(b), apex(a), radius(r) {}

        T height() const { return distance(base, apex); }
    };

    /// AABB (軸平行境界ボックス)
    template<typename T>
    struct AABB {
        Point3D<T> min, max;

        AABB() = default;
        AABB(const Point3D<T>& mn, const Point3D<T>& mx) : min(mn), max(mx) {}

        bool contains(const Point3D<T>& p) const {
            return p.x >= min.x && p.x <= max.x
                && p.y >= min.y && p.y <= max.y
                && p.z >= min.z && p.z <= max.z;
        }

        bool intersects(const AABB& other) const {
            return min.x <= other.max.x && max.x >= other.min.x
                && min.y <= other.max.y && max.y >= other.min.y
                && min.z <= other.max.z && max.z >= other.min.z;
        }

        Point3D<T> center() const {
            return {(min.x + max.x) / T{2}, (min.y + max.y) / T{2}, (min.z + max.z) / T{2}};
        }
    };

    /// トーラス (Y軸回りの回転体)
    template<typename T>
    struct Torus {
        Point3D<T> center;
        T majorRadius{};  // 中心円の半径 R
        T minorRadius{};  // 管の半径 r

        Torus() = default;
        Torus(const Point3D<T>& c, T R, T r)
            : center(c), majorRadius(R), minorRadius(r) {}
    };

    /// ディスク (円盤 — 法線 + 中心 + 内外半径)
    template<typename T>
    struct Disc {
        Point3D<T> center;
        Point3D<T> normal;
        T outerRadius{};
        T innerRadius{};  // 0 = 穴なし

        Disc() = default;
        Disc(const Point3D<T>& c, const Point3D<T>& n, T outer, T inner = T(0))
            : center(c), normal(n.normalized()), outerRadius(outer), innerRadius(inner) {}
    };

    // ストリーム出力
    template<typename T>
    std::ostream& operator<<(std::ostream& os, const Point2D<T>& p) {
        return os << "(" << p.x << ", " << p.y << ")";
    }

    template<typename T>
    std::ostream& operator<<(std::ostream& os, const Point3D<T>& p) {
        return os << "(" << p.x << ", " << p.y << ", " << p.z << ")";
    }

} // namespace calx::image

#endif // CALX_IMAGE_GEOMETRY_PRIMITIVES_HPP
