// Copyright (c) 2005 Rijksuniversiteit Groningen (Netherlands)
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
// You can redistribute it and/or modify it under the terms of the GNU
// General Public License as published by the Free Software Foundation,
// either version 3 of the License, or (at your option) any later version.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0+
//
//
// Author(s)     : Nico Kruithof <Nico@cs.rug.nl>

#ifndef CGAL_SKIN_SURFACE_QUADRATIC_SURFACE_H
#define CGAL_SKIN_SURFACE_QUADRATIC_SURFACE_H

#include <CGAL/license/Skin_surface_3.h>

#include <CGAL/Skin_surface_traits_3.h>

namespace CGAL {

template < class SkinSurfaceQuadraticSurfaceTraits_3 >
class Skin_surface_quadratic_surface_3
{
public:
  typedef SkinSurfaceQuadraticSurfaceTraits_3     K;
  typedef typename K::FT                          FT;
  typedef typename K::Point_3                     Point;
  typedef typename K::Vector_3                    Vector;
  typedef typename K::Segment_3                   Segment;
  typedef typename K::Weighted_point_3            Weighted_point;

  Skin_surface_quadratic_surface_3()
    : dim(-1), p(0,0,0), c(0)
  {
    for (int i=0; i<6; i++) Q[i] = 0;
  }

  Skin_surface_quadratic_surface_3(FT Qinput[], Point p, FT c, int d)
    : dim(10+d), p(p), c(c)
  {
    for (int i=0; i<6; i++) Q[i] = Qinput[i];
  }

  Skin_surface_quadratic_surface_3(Weighted_point wp0, FT s)
    : dim(0), p(wp0.point()), c(-s*(1-s)*wp0.weight())
  {
    CGAL_PROFILER(std::string("Constructor : ") +
                  std::string(CGAL_PRETTY_FUNCTION));
    Q[1] = Q[3] = Q[4] = 0;
    Q[0] = Q[2] = Q[5] = (1-s);
  }

  Skin_surface_quadratic_surface_3(Weighted_point wp0,
                                   Weighted_point wp1,
                                   FT s) : dim(1)
  {
    CGAL_PROFILER(std::string("Constructor : ") +
                  std::string(CGAL_PRETTY_FUNCTION));

    K k;
    p = k.construct_weighted_circumcenter_3_object()(wp0,wp1);
    c = s*(1-s)*k.compute_squared_radius_smallest_orthogonal_sphere_3_object()(wp0,wp1);
    Vector t = k.construct_vector_3_object()(k.construct_point_3_object()(wp1),
                                             k.construct_point_3_object()(wp0));

    FT den = t*t;
    Q[0] = (-  t.x()*t.x()/den + (1-s));

    Q[1] = (-2*t.y()*t.x()/den);
    Q[2] = (-  t.y()*t.y()/den + (1-s));

    Q[3] = (-2*t.z()*t.x()/den);
    Q[4] = (-2*t.z()*t.y()/den);
    Q[5] = (-  t.z()*t.z()/den + (1-s));
  }

  Skin_surface_quadratic_surface_3(Weighted_point wp0,
                                   Weighted_point wp1,
                                   Weighted_point wp2,
                                   FT s) : dim(2)
  {
    CGAL_PROFILER(std::string("Constructor : ") +
                  std::string(CGAL_PRETTY_FUNCTION));

    K k;
    p = k.construct_weighted_circumcenter_3_object()(wp0,wp1,wp2);
    c = s*(1-s)*k.compute_squared_radius_smallest_orthogonal_sphere_3_object()(wp0,wp1,wp2);

    Vector t = K().construct_orthogonal_vector_3_object()(k.construct_point_3_object()(wp0),
                                                          k.construct_point_3_object()(wp1),
                                                          k.construct_point_3_object()(wp2));

    FT den = t*t;
    Q[0] = -(-  t.x()*t.x()/den + s);

    Q[1] = -(-2*t.y()*t.x()/den);
    Q[2] = -(-  t.y()*t.y()/den + s);

    Q[3] = -(-2*t.z()*t.x()/den);
    Q[4] = -(-2*t.z()*t.y()/den);
    Q[5] = -(-  t.z()*t.z()/den + s);

  }
  Skin_surface_quadratic_surface_3(Weighted_point wp0,
                                   Weighted_point wp1,
                                   Weighted_point wp2,
                                   Weighted_point wp3,
                                   FT s) : dim(3)
  {
    CGAL_PROFILER(std::string("Constructor : ") +
                  std::string(CGAL_PRETTY_FUNCTION));

    K k;
    p = k.construct_weighted_circumcenter_3_object()(wp0,wp1,wp2,wp3);
    c = s*(1-s)*k.compute_squared_radius_smallest_orthogonal_sphere_3_object()(wp0,wp1,wp2,wp3);
    Q[1] = Q[3] = Q[4] = 0;
    Q[0] = Q[2] = Q[5] = -s;
  }

  template <class Input_point>
  FT value(Input_point const &x) const
  {
    typedef Cartesian_converter<typename Input_point::R, K> Converter;

    FT vx = Converter()(x.x()) - p.x();
    FT vy = Converter()(x.y()) - p.y();
    FT vz = Converter()(x.z()) - p.z();

    return vx*(Q[0]*vx) +
           vy*(Q[1]*vx+Q[2]*vy) +
           vz*(Q[3]*vx+Q[4]*vy+Q[5]*vz) +
           c;
  }
  template <class Input_point>
  Sign sign(Input_point const &x) const {
    return CGAL_NTS sign(value(x));
  }

  template <class Input_point>
  typename Input_point::R::Vector_3 gradient(Input_point const &x)
  {
    // NGHK: We use the kernel trick again: DOCUMENT!!!!
    typedef Cartesian_converter<typename Input_point::R, K> Converter;
    typedef Cartesian_converter<K, typename Input_point::R> Inv_converter;
    Converter conv;
    Inv_converter inv_conv;
    Point xp = conv(x);
    return inv_conv(compute_gradient(xp));
  }

  /// Construct the intersection point with the segment (p0,p1)
  template <class Input_point>
  Input_point to_surface(const Input_point &p0, const Input_point &p1)
  {
    // NGHK: We use the kernel trick again: DOCUMENT!!!!
    typedef Cartesian_converter<typename Input_point::R, K> Converter;
    Converter conv;

    Input_point pp0 = p0, pp1 = p1, mid;
    double sq_d = to_double(squared_distance(pp0,pp1));

    if (value(conv(pp1)) < value(conv(pp0))) {
      std::swap(pp0, pp1);
    }

    while (sq_d > 1.e-10) {
      mid = midpoint(pp0,pp1);
      if (value(conv(mid)) > 0) {
        pp1 = mid;
      } else {
        pp0 = mid;
      }
      sq_d /= 4;
    }
    return midpoint(pp0,pp1);
  }

private:
  //template <class Input_point>
  Vector compute_gradient(Point const &x)
  {
    FT vx = x.x() - p.x();
    FT vy = x.y() - p.y();
    FT vz = x.z() - p.z();

    return Vector(2*Q[0]*vx + Q[1]*vy + Q[3]*vz,
                  Q[1]*vx + 2*Q[2]*vy + Q[4]*vz,
                  Q[3]*vx + Q[4]*vy + 2*Q[5]*vz);
  }

  int dim;
  FT Q[6];
  Point p;
  FT c;
};

} //namespace CGAL

#endif // CGAL_SKIN_SURFACE_QUADRATIC_SURFACE_H