// Copyright (c) 2008 INRIA(France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org)
//
// $URL: https://github.com/CGAL/cgal/blob/v6.2-beta1/Intersections_3/include/CGAL/Intersections_3/internal/Plane_3_Plane_3_Plane_3_intersection.h $
// $Id: include/CGAL/Intersections_3/internal/Plane_3_Plane_3_Plane_3_intersection.h 431cfb5357d $
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
//
//
// Author(s)     : Geert-Jan Giezeman,
//                 Sylvain Pion

#ifndef CGAL_INTERNAL_INTERSECTIONS_3_PLANE_3_PLANE_3_PLANE_3_INTERSECTION_H
#define CGAL_INTERNAL_INTERSECTIONS_3_PLANE_3_PLANE_3_PLANE_3_INTERSECTION_H

#include <CGAL/Intersection_traits_3.h>
#include <CGAL/Intersections_3/internal/Plane_3_Plane_3_intersection.h>
#include <CGAL/Intersections_3/internal/Line_3_Plane_3_intersection.h>

#include <optional>
#include <variant>

namespace CGAL {
namespace Intersections {
namespace internal {

//  triple plane intersection
template <class K>
std::optional<typename K::Point_3>
intersection_point_RT(const typename K::Plane_3& plane1,
                      const typename K::Plane_3& plane2,
                      const typename K::Plane_3& plane3,
                      const K&)
{
  typedef typename K::RT RT;

  const RT& m00 = plane1.a();
  const RT& m01 = plane1.b();
  const RT& m02 = plane1.c();
  const RT& b0  = - plane1.d();
  const RT& m10 = plane2.a();
  const RT& m11 = plane2.b();
  const RT& m12 = plane2.c();
  const RT& b1  = - plane2.d();
  const RT& m20 = plane3.a();
  const RT& m21 = plane3.b();
  const RT& m22 = plane3.c();
  const RT& b2  = - plane3.d();

  // Minors common to two determinants
  const RT minor_0 = m00*m11 - m10*m01;
  const RT minor_1 = m00*m21 - m20*m01;
  const RT minor_2 = m10*m21 - m20*m11;

  const RT den = minor_0*m22 - minor_1*m12 + minor_2*m02; // determinant of M

  if(is_zero(den)){
    return std::nullopt;
  }

  const RT num3 = minor_0*b2 - minor_1*b1 + minor_2*b0;  // determinant of M with M[x:2] swapped with [b0,b1,b2]

  // Minors common to two determinants
  const RT minor_3 = b0*m12 - b1*m02;
  const RT minor_4 = b0*m22 - b2*m02;
  const RT minor_5 = b1*m22 - b2*m12;

  // num1 has opposite signs because b0 and M[:1] have been swapped
  const RT num1 = - minor_3*m21 + minor_4*m11 - minor_5*m01;  // determinant of M with M[x:0] swapped with [b0,b1,b2]
  const RT num2 = minor_3*m20 - minor_4*m10 + minor_5*m00;  // determinant of M with M[x:1] swapped with [b0,b1,b2]

  return typename K::Point_3(num1, num2, num3, den);
}

template <class K>
std::optional<typename K::Point_3>
intersection_point(const typename K::Plane_3& plane1,
                   const typename K::Plane_3& plane2,
                   const typename K::Plane_3& plane3,
                   const K& k)
{
  return intersection_point_RT(plane1, plane2, plane3, k);
}

template <class K>
std::optional<std::variant<typename K::Point_3, typename K::Line_3, typename K::Plane_3> >
intersection(const typename K::Plane_3& plane1,
             const typename K::Plane_3& plane2,
             const typename K::Plane_3& plane3,
             const K& k)
{
  typedef typename std::optional<std::variant<typename K::Point_3,
                                                  typename K::Line_3,
                                                  typename K::Plane_3> > result_type;

  typedef typename K::Line_3       Line_3;
  typedef typename K::Plane_3      Plane_3;

  auto res = intersection_point(plane1,plane2,plane3, k);
  if (res)
    return result_type(*res);

  // Intersection between plane1 and plane2 can either be
  // a line, a plane, or empty.
  typename Intersection_traits<K, Plane_3, Plane_3>::result_type
      o12 = internal::intersection(plane1, plane2, k);

  if(o12)
  {
    if(const Line_3* l = intersect_get<Line_3>(o12))
    {
      if (internal::do_intersect(*l, plane3, k))
        return result_type(*l);
    }
    else
    {
      CGAL_assertion(intersect_get<Plane_3>(o12) != nullptr);
      // either line or plane
      typename Intersection_traits<K, Plane_3, Plane_3>::result_type
          v = internal::intersection(plane3, plane1, k);
      if(v)
      {
        if( intersect_get<Plane_3>(v)!=nullptr)
          return result_type(plane1);
        else if(const Line_3* l = intersect_get<Line_3>(v))
          return result_type(*l);
      }
    }
  }

  return result_type();
}

} // namespace internal
} // namespace Intersections
} // namespace CGAL

#endif // CGAL_INTERNAL_INTERSECTIONS_3_PLANE_3_PLANE_3_PLANE_3_INTERSECTION_H
