// Copyright (c) 2008  INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL: https://github.com/CGAL/cgal/blob/v6.1.2/Circular_kernel_3/include/CGAL/Circular_kernel_3/internal_functions_on_sphere_3.h $
// $Id: include/CGAL/Circular_kernel_3/internal_functions_on_sphere_3.h 37e2b278910 $
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Monique Teillaud, Sylvain Pion, Pedro Machado

// Partially supported by the IST Programme of the EU as a
// STREP (FET Open) Project under Contract No  IST-006413
// (ACS -- Algorithms for Complex Shapes)

#ifndef CGAL_SPHERICAL_KERNEL_PREDICATES_ON_SPHERE_3_H
#define CGAL_SPHERICAL_KERNEL_PREDICATES_ON_SPHERE_3_H

#include <CGAL/license/Circular_kernel_3.h>


#include <CGAL/Circular_kernel_3/Intersection_traits.h>
#include <CGAL/number_utils.h>
#include <utility>
#include <vector>

namespace CGAL {
  namespace SphericalFunctors {

    template < class SK >
    typename SK::Algebraic_kernel::Polynomials_for_line_3
    get_equation( const typename SK::Line_3 & l)
    {
      typedef typename SK::Algebraic_kernel Algebraic_kernel;
      return Algebraic_kernel().construct_polynomials_for_line_3_object()
          (l.to_vector().x(), l.point().x(),
           l.to_vector().y(), l.point().y(),
           l.to_vector().z(), l.point().z());
    }

    template < class SK >
    typename SK::Polynomial_1_3
    get_equation( const typename SK::Plane_3 & s )
    {
      typedef typename SK::Algebraic_kernel Algebraic_kernel;
      return Algebraic_kernel().construct_polynomial_1_3_object()
          ( s.a(), s.b(), s.c(), s.d() );
    }

    template < class SK >
    typename SK::Polynomial_for_spheres_2_3
    get_equation( const typename SK::Sphere_3 & s )
    {
      typedef typename SK::Point_3    Point_3;
      typedef typename SK::Algebraic_kernel   Algebraic_kernel;
      Point_3 center = s.center();
      return Algebraic_kernel().construct_polynomial_for_spheres_2_3_object()
          ( center.x(), center.y(), center.z(), s.squared_radius() );
    }

    template < class SK >
    typename SK::Polynomials_for_circle_3
    get_equation( const typename SK::Circle_3 & c )
    {
      typedef typename SK::Algebraic_kernel Algebraic_kernel;
      return std::make_pair ( Algebraic_kernel().construct_polynomial_for_spheres_2_3_object()
                               (c.center().x(), c.center().y(),
                                c.center().z(), c.squared_radius()),
                              Algebraic_kernel().construct_polynomial_1_3_object()
                               (c.supporting_plane().a(),
                                c.supporting_plane().b(),
                                c.supporting_plane().c(),
                                c.supporting_plane().d()));
    }

    template < class SK >
    typename SK::Sphere_3
    construct_sphere_3(const typename SK::Polynomial_for_spheres_2_3 &eq)
    {
      typedef typename SK::Sphere_3 Sphere_3;
      typedef typename SK::Point_3  Point_3;
      return Sphere_3(Point_3(eq.a(),eq.b(),eq.c()),eq.r_sq());
    }


    template < class SK >
    inline
    typename SK::Bounded_side
    bounded_side(const typename SK::Sphere_3 &s,
                 const typename SK::Circular_arc_point_3 &p) {
      typedef typename SK::Algebraic_kernel Algebraic_kernel;
      typedef typename SK::Polynomial_for_spheres_2_3 Equation;
      Equation equation = get_equation<SK>(s);
      Sign sign = Algebraic_kernel().sign_at_object()(equation,p.rep().coordinates());
      if(sign == NEGATIVE) return ON_BOUNDED_SIDE;
      else if(sign == POSITIVE) return ON_UNBOUNDED_SIDE;
      else return ON_BOUNDARY;
    }

    template < class SK >
    inline
    typename SK::Bounded_side
    bounded_side(const typename SK::Circle_3 &c,
                 const typename SK::Circular_arc_point_3 &p) {
      typedef typename SK::Algebraic_kernel Algebraic_kernel;
      typedef typename SK::Polynomial_for_spheres_2_3 Equation;
      CGAL_kernel_assertion(SK().has_on_3_object()(c.supporting_plane(),p));
      Equation equation = get_equation<SK>(c.diametral_sphere());
      Sign sign = Algebraic_kernel().sign_at_object()(equation,p.rep().coordinates());
      if(sign == NEGATIVE) return ON_BOUNDED_SIDE;
      else if(sign == POSITIVE) return ON_UNBOUNDED_SIDE;
      else return ON_BOUNDARY;
    }

    template < class SK >
    inline
    bool
    non_oriented_equal(const typename SK::Sphere_3 & s1,
                       const typename SK::Sphere_3 & s2)
    {
      // Should we compare anyway even if they are degenerated?
      CGAL_kernel_assertion(!(s1.is_degenerate() || s2.is_degenerate()));
      return s1.center() == s2.center() &&
      s1.squared_radius() == s2.squared_radius();
    }

    template < class SK >
    inline
    bool
    non_oriented_equal(const typename SK::Plane_3 & p1,
                       const typename SK::Plane_3 & p2)
    {
      // Should we compare anyway even if they are degenerated?
      CGAL_kernel_assertion(!(p1.is_degenerate() || p2.is_degenerate()));
      if(CGAL::is_zero(p1.a())) {
        if(!CGAL::is_zero(p2.a())) return false;
        if(CGAL::is_zero(p1.b())) {
          if(!CGAL::is_zero(p2.b())) return false;
          return p1.c() * p2.d() == p1.d() * p2.c();
        }
        return (p2.c() * p1.b() == p1.c() * p2.b()) &&
               (p2.d() * p1.b() == p1.d() * p2.b());
      }
      return (p2.b() * p1.a() == p1.b() * p2.a()) &&
             (p2.c() * p1.a() == p1.c() * p2.a()) &&
             (p2.d() * p1.a() == p1.d() * p2.a());
    }

    template < class SK >
    inline
    bool
    non_oriented_equal(const typename SK::Circle_3 & c1,
                       const typename SK::Circle_3 & c2) {
      // We see degeneracies on the other non_oriented_equal functions
      if(!non_oriented_equal<SK>(c1.diametral_sphere(), c2.diametral_sphere())) return false;
      if(!non_oriented_equal<SK>(c1.supporting_plane(), c2.supporting_plane())) return false;
      return true;
    }

    template< class SK>
    bool
    non_oriented_equal(const typename SK::Line_3 &l1,
                       const typename SK::Line_3 &l2)
    {
      typedef typename SK::Vector_3 Vector_3;
      if(!SK().has_on_3_object()(l1, l2.point())) return false;

      const Vector_3& v1 = l1.to_vector();
      const Vector_3& v2 = l2.to_vector();

      if(v1.x() * v2.y() != v1.y() * v2.x()) return false;
      if(v1.x() * v2.z() != v1.z() * v2.x()) return false;
      if(v1.y() * v2.z() != v1.z() * v2.y()) return false;

      return true;
    }

    template< class SK>
    bool
    non_oriented_equal( const typename SK::Line_arc_3 &l1,
                        const typename SK::Line_arc_3 &l2)
    {
      if(!non_oriented_equal<SK>(l1.supporting_line(),
                                 l2.supporting_line())) return false;
      return (l1.lower_xyz_extremity() == l2.lower_xyz_extremity()) &&
             (l1.higher_xyz_extremity() == l2.higher_xyz_extremity());
    }

    template< class SK>
    bool
    non_oriented_equal( const typename SK::Circular_arc_3 &c1,
                        const typename SK::Circular_arc_3 &c2)
    {
      if(!non_oriented_equal<SK>(c1.supporting_circle(),
                                 c2.supporting_circle())) return false;
      if(c1.rep().is_full() && c2.rep().is_full())
        return true;
      return (c1.source() == c2.source()) &&
             (c1.target() == c2.target());
    }

   namespace internal {
     // we need to pass the result_type, to hack around the fact that
     // object doesn't support operator=(const T&) and so we keep backwards compatibility
     template<typename SK>
     struct pair_transform {
       auto operator()(const std::pair< typename SK::Root_for_spheres_2_3, unsigned >& p) {
         return std::make_pair(typename SK::Circular_arc_point_3(p.first), p.second);
       }
     };

   // obscure trick: calculate the regular intersection between two
   // objects, throw it away if empty, else check if the result was a
   // point_3, if so turn in into a (circular_arc_point, unsigned)
   // pair, otherwise just dump the value with conversion to RT
   // (again: converting to RT before assigning to the Iterator is
   // just to keep object working)
   template<typename SK, typename OutputIterator>
   struct Point_conversion_visitor {
     Point_conversion_visitor(const OutputIterator& it) : it(it) {}
     template<typename T>
     OutputIterator operator()(const T& t) { *it++ = t; return it; }

     OutputIterator operator()(const typename SK::Point_3& p) {
       // 2 multiplicities
       *it++ = std::make_pair(typename SK::Circular_arc_point_3(p), 2u);
       return it;
     }
     OutputIterator it;
   };

   } // namespace internal

    template < class SK, class OutputIterator >
    OutputIterator
    intersect_3(const typename SK::Sphere_3 & s,
                const typename SK::Line_3 & l,
                OutputIterator res)
    {
      typedef typename SK::Algebraic_kernel                          Algebraic_kernel;
      typedef typename SK::Polynomial_for_spheres_2_3  Equation_sphere;
      typedef typename SK::Polynomials_for_line_3      Equation_line;
      typedef typename SK::Root_for_spheres_2_3        Root_for_spheres_2_3;
      CGAL_kernel_precondition(!s.is_degenerate());
      CGAL_kernel_precondition(!l.is_degenerate());
      Equation_sphere e1 = get_equation<SK>(s);
      Equation_line e2 = get_equation<SK>(l);
      typedef std::vector< std::pair < Root_for_spheres_2_3, unsigned > >
        solutions_container;
      solutions_container solutions;
      Algebraic_kernel().solve_object()(e1, e2, std::back_inserter(solutions));

      return std::transform(solutions.begin(), solutions.end(), res, internal::pair_transform<SK>());
    }

    // The special 3 object functions
    template < class SK, class OutputIterator >
    OutputIterator
    intersect_3(const typename SK::Sphere_3 & s1,
                const typename SK::Sphere_3 & s2,
                const typename SK::Sphere_3 & s3,
                OutputIterator res)
    {
       // circle_3, sphere_3, (circular_arc_point_3, uint)
       typedef typename SK::Polynomial_for_spheres_2_3  Equation_sphere;
       typedef typename SK::Root_for_spheres_2_3  Root_for_spheres_2_3;
       typedef typename SK::Circular_arc_point_3  Circular_arc_point_3;
       typedef typename SK::Circle_3  Circle_3;
       typedef typename SK::Point_3  Point_3;
       typedef typename SK::Sphere_3  Sphere_3;
       typedef typename SK::Algebraic_kernel  Algebraic_kernel;
       typedef typename SK::Boolean Bool;

       CGAL_kernel_precondition(!s1.is_degenerate());
       CGAL_kernel_precondition(!s2.is_degenerate());
       CGAL_kernel_precondition(!s3.is_degenerate());
       if(non_oriented_equal<SK>(s1,s2) && non_oriented_equal<SK>(s2,s3)) {
         *res++ = s1;
         return res;
       }
       if(non_oriented_equal<SK>(s1,s2)) {
         if(typename Intersection_traits<SK, Sphere_3, Sphere_3>::result_type v =
            SK().intersect_3_object()(s1, s3)) {
           internal::Point_conversion_visitor<SK, OutputIterator> visitor(res);
           return std::visit(visitor,
             *v);
         }
         return res;
       }
       if(non_oriented_equal<SK>(s1,s3) || non_oriented_equal<SK>(s2,s3)) {
         if(typename Intersection_traits<SK, Sphere_3, Sphere_3>::result_type v =
            SK().intersect_3_object()(s1, s2)) {
           internal::Point_conversion_visitor<SK, OutputIterator> visitor(res);
           return std::visit(
             visitor,
             *v);
         }
         return res;
       }
       if(SK().collinear_3_object()(s1.center(),s2.center(),s3.center())) {
         typename Intersection_traits<SK, Sphere_3, Sphere_3>::result_type v =
           SK().intersect_3_object()(s1, s2);
         if(!v) return res;
         if(const Point_3* p = CGAL::Intersections::internal::intersect_get<Point_3>(v)) {
            if(SK().has_on_3_object()(s3, *p)) {
              *res++ = std::make_pair(Circular_arc_point_3(*p),2u);
            }
             return res;
         }
         if(const Circle_3* c = CGAL::Intersections::internal::intersect_get<Circle_3>(v)) {
            if(static_cast<Bool>(SK().has_on_3_object()(s3, *c))) {
              *res++ = *c;
            }
           return res;
         }
         return res;
       }
       Equation_sphere e1 = get_equation<SK>(s1);
       Equation_sphere e2 = get_equation<SK>(s2);
       Equation_sphere e3 = get_equation<SK>(s3);
       typedef std::vector< std::pair < Root_for_spheres_2_3, unsigned > >
         algebraic_solutions_container;
       algebraic_solutions_container solutions;
       Algebraic_kernel().solve_object()(e1, e2, e3, std::back_inserter(solutions));

       return std::transform(solutions.begin(), solutions.end(), res, internal::pair_transform<SK>());
    }

    template < class SK, class OutputIterator >
    OutputIterator
    intersect_3(const typename SK::Plane_3 & p,
                const typename SK::Sphere_3 & s1,
                const typename SK::Sphere_3 & s2,
                OutputIterator res)
    {
      typedef typename SK::Root_for_spheres_2_3  Root_for_spheres_2_3;
      typedef typename SK::Polynomial_for_spheres_2_3  Equation_sphere;
      typedef typename SK::Polynomial_1_3  Equation_plane;
      typedef typename SK::Plane_3  Plane_3;
      typedef typename SK::Sphere_3 Sphere_3;
      typedef typename SK::Algebraic_kernel  Algebraic_kernel;
      CGAL_kernel_precondition(!p.is_degenerate());
      CGAL_kernel_precondition(!s1.is_degenerate());
      CGAL_kernel_precondition(!s2.is_degenerate());
      if(non_oriented_equal<SK>(s1,s2)) {
        if(typename Intersection_traits<SK, Plane_3, Sphere_3>::result_type v =
            SK().intersect_3_object()(p, s1)) {
           internal::Point_conversion_visitor<SK, OutputIterator> visitor(res);
           return std::visit(
             visitor,
             *v);
         }
         return res;
      }
      Plane_3 radical_p = SK().construct_radical_plane_3_object()(s1,s2);
      if(non_oriented_equal<SK>(p,radical_p)) {
        if(typename Intersection_traits<SK, Plane_3, Sphere_3>::result_type v =
            SK().intersect_3_object()(p, s1)) {
           internal::Point_conversion_visitor<SK, OutputIterator> visitor(res);
           return std::visit(
             visitor,
             *v);
         }
         return res;
      }
      Equation_sphere e1 = get_equation<SK>(s1);
      Equation_sphere e2 = get_equation<SK>(s2);
      Equation_plane e3 = get_equation<SK>(p);
      typedef std::vector< std::pair < Root_for_spheres_2_3, unsigned > >
        algebraic_solutions_container;
      algebraic_solutions_container solutions;
      Algebraic_kernel().solve_object()(e1, e2, e3, std::back_inserter(solutions));
      return std::transform(solutions.begin(), solutions.end(), res, internal::pair_transform<SK>());
    }

    template < class SK, class OutputIterator >
    OutputIterator
    intersect_3(const typename SK::Plane_3 & p1,
                const typename SK::Plane_3 & p2,
                const typename SK::Sphere_3 & s,
                OutputIterator res)
    {
      typedef typename SK::Root_for_spheres_2_3  Root_for_spheres_2_3;
      typedef typename SK::Polynomial_for_spheres_2_3  Equation_sphere;
      typedef typename SK::Polynomial_1_3  Equation_plane;
      typedef typename SK::Algebraic_kernel  Algebraic_kernel;
      typedef typename SK::Plane_3 Plane_3;
      typedef typename SK::Sphere_3 Sphere_3;
      CGAL_kernel_precondition(!p1.is_degenerate());
      CGAL_kernel_precondition(!p2.is_degenerate());
      CGAL_kernel_precondition(!s.is_degenerate());
      if(non_oriented_equal<SK>(p1,p2)) {
        if(typename Intersection_traits<SK, Plane_3, Sphere_3>::result_type v =
            SK().intersect_3_object()(p1, s)) {
           internal::Point_conversion_visitor<SK, OutputIterator> visitor(res);
           return std::visit(
             visitor,
             *v);
         }
         return res;
      }
      Equation_plane e1 = get_equation<SK>(p1);
      Equation_plane e2 = get_equation<SK>(p2);
      Equation_sphere e3 = get_equation<SK>(s);
      typedef std::vector< std::pair < Root_for_spheres_2_3, unsigned > >
        algebraic_solutions_container;
      algebraic_solutions_container solutions;
      Algebraic_kernel().solve_object()(e1, e2, e3, std::back_inserter(solutions));
      return std::transform(solutions.begin(), solutions.end(), res, internal::pair_transform<SK>());
    }

    template < class SK, class OutputIterator >
    OutputIterator
    intersect_3(const typename SK::Circle_3 & c,
                const typename SK::Plane_3 & p,
                OutputIterator res)
    {
      return intersect_3<SK>(p,c.supporting_plane(),c.diametral_sphere(),res);
    }

    template < class SK, class OutputIterator >
    OutputIterator
    intersect_3(const typename SK::Circle_3 & c,
                const typename SK::Sphere_3 & s,
                OutputIterator res)
    {
      return intersect_3<SK>(c.supporting_plane(),s,c.diametral_sphere(),res);
    }

    template < class SK, class OutputIterator >
    OutputIterator
    intersect_3(const typename SK::Circle_3 & c1,
                const typename SK::Circle_3 & c2,
                OutputIterator res)
    {
      typedef typename SK::Root_for_spheres_2_3     Root_for_spheres_2_3;
      typedef typename SK::Polynomials_for_circle_3 Equation_circle;
      typedef typename SK::Algebraic_kernel         Algebraic_kernel;

      if(non_oriented_equal<SK>(c1,c2)) {
         *res++ = c1;
         return res;
      }
      Equation_circle e1 = get_equation<SK>(c1);
      Equation_circle e2 = get_equation<SK>(c2);
      typedef std::vector< std::pair < Root_for_spheres_2_3, unsigned > >
        algebraic_solutions_container;
      algebraic_solutions_container solutions;
      Algebraic_kernel().solve_object()(e1, e2, std::back_inserter(solutions));
      return std::transform(solutions.begin(), solutions.end(), res, internal::pair_transform<SK>());
    }

    template < class SK, class OutputIterator >
    OutputIterator
    intersect_3(const typename SK::Circle_3 & c,
                const typename SK::Line_3 & l,
                OutputIterator res)
    {
      typedef typename SK::Root_for_spheres_2_3     Root_for_spheres_2_3;
      typedef typename SK::Polynomials_for_circle_3 Equation_circle;
      typedef typename SK::Polynomials_for_line_3   Equation_line;

      typedef typename SK::Algebraic_kernel  Algebraic_kernel;
      CGAL_kernel_precondition(!l.is_degenerate());
      Equation_circle e1 = get_equation<SK>(c);
      Equation_line e2 = get_equation<SK>(l);
      typedef std::vector< std::pair < Root_for_spheres_2_3, unsigned > >
        algebraic_solutions_container;
      algebraic_solutions_container solutions;
      Algebraic_kernel().solve_object()(e1, e2, std::back_inserter(solutions));
      return std::transform(solutions.begin(), solutions.end(), res, internal::pair_transform<SK>());
    }

    // At the moment we don't need those functions
    // But in the future maybe (some make_x_monotone? etc..)
    template <class SK>
    typename SK::Circular_arc_point_3
    x_extremal_point(const typename SK::Sphere_3 & c, bool i)
    {
      typedef typename SK::Algebraic_kernel   Algebraic_kernel;
      return Algebraic_kernel().x_critical_points_object()(typename SK::Get_equation()(c),i);
    }

    template <class SK,class OutputIterator>
    OutputIterator
    x_extremal_points(const typename SK::Sphere_3 & c, OutputIterator res)
    {
      typedef typename SK::Algebraic_kernel   Algebraic_kernel;
      return Algebraic_kernel().x_critical_points_object()(typename SK::Get_equation()(c),res);
    }

    template <class SK>
    typename SK::Circular_arc_point_3
    y_extremal_point(const typename SK::Sphere_3 & c, bool i)
    {
      typedef typename SK::Algebraic_kernel   Algebraic_kernel;
      return Algebraic_kernel().y_critical_points_object()(typename SK::Get_equation()(c),i);
    }

    template <class SK,class OutputIterator>
    OutputIterator
    y_extremal_points(const typename SK::Sphere_3 & c, OutputIterator res)
    {
      typedef typename SK::Algebraic_kernel   Algebraic_kernel;
      return Algebraic_kernel().y_critical_points_object()(typename SK::Get_equation()(c),res);
    }

    template <class SK>
    typename SK::Circular_arc_point_3
    z_extremal_point(const typename SK::Sphere_3 & c, bool i)
    {
      typedef typename SK::Algebraic_kernel   Algebraic_kernel;
      return Algebraic_kernel().z_critical_points_object()(typename SK::Get_equation()(c),i);
    }

    template <class SK,class OutputIterator>
    OutputIterator
    z_extremal_points(const typename SK::Sphere_3 & c, OutputIterator res)
    {
      typedef typename SK::Algebraic_kernel   Algebraic_kernel;
      return Algebraic_kernel().z_critical_points_object()(typename SK::Get_equation()(c),res);
    }

    template <class SK>
    typename SK::Circular_arc_point_3
    x_extremal_point(const typename SK::Circle_3 & c, bool i)
    {
      typedef typename SK::Algebraic_kernel   Algebraic_kernel;
      return Algebraic_kernel().x_critical_points_object()(typename SK::Get_equation()(c),i);
    }

    template <class SK,class OutputIterator>
    OutputIterator
    x_extremal_points(const typename SK::Circle_3 & c, OutputIterator res)
    {
      typedef typename SK::Algebraic_kernel   Algebraic_kernel;
      return Algebraic_kernel().x_critical_points_object()(typename SK::Get_equation()(c),res);
    }

    template <class SK>
    typename SK::Circular_arc_point_3
    y_extremal_point(const typename SK::Circle_3 & c, bool i)
    {
      typedef typename SK::Algebraic_kernel   Algebraic_kernel;
      return Algebraic_kernel().y_critical_points_object()(typename SK::Get_equation()(c),i);
    }

    template <class SK,class OutputIterator>
    OutputIterator
    y_extremal_points(const typename SK::Circle_3 & c, OutputIterator res)
    {
      typedef typename SK::Algebraic_kernel   Algebraic_kernel;
      return Algebraic_kernel().y_critical_points_object()(typename SK::Get_equation()(c),res);
    }

    template <class SK>
    typename SK::Circular_arc_point_3
    z_extremal_point(const typename SK::Circle_3 & c, bool i)
    {
      typedef typename SK::Algebraic_kernel   Algebraic_kernel;
      return Algebraic_kernel().z_critical_points_object()(typename SK::Get_equation()(c),i);
    }

    template <class SK,class OutputIterator>
    OutputIterator
    z_extremal_points(const typename SK::Circle_3 & c, OutputIterator res)
    {
      typedef typename SK::Algebraic_kernel   Algebraic_kernel;
      return Algebraic_kernel().z_critical_points_object()(typename SK::Get_equation()(c),res);
    }

  }//SphericalFunctors
}//CGAL

#endif //CGAL_SPHERICAL_KERNEL_PREDICATES_ON_SPHERE_3_H
