// 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$
// 
//
// Author(s)     : Nico Kruithof <Nico@cs.rug.nl>

#ifndef CGAL_SKIN_SURFACE_QUADRATIC_SURFACE_H
#define CGAL_SKIN_SURFACE_QUADRATIC_SURFACE_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 CGAL::Weighted_point<Point,FT>  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];
  }
  int dim;
  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));
    Regular_triangulation_euclidean_traits_3<K> reg_traits;
    p = reg_traits.construct_weighted_circumcenter_3_object()(wp0,wp1);
    c = s*(1-s)*reg_traits.compute_squared_radius_smallest_orthogonal_sphere_3_object()(wp0,wp1);
    Vector t = wp0-wp1;

    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));
    Regular_triangulation_euclidean_traits_3<K> reg_traits;
    p = reg_traits.construct_weighted_circumcenter_3_object()(wp0,wp1,wp2);
    c = s*(1-s)*reg_traits.compute_squared_radius_smallest_orthogonal_sphere_3_object()(wp0,wp1,wp2);
    
    Vector t = K().construct_orthogonal_vector_3_object()(wp0,wp1,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));
    Regular_triangulation_euclidean_traits_3<K> reg_traits;
    p = reg_traits.construct_weighted_circumcenter_3_object()(wp0,wp1,wp2,wp3);
    c = s*(1-s)*reg_traits.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);
  }
	
  FT Q[6]; 
  Point p; 
  FT c;
};

} //namespace CGAL

#endif // CGAL_SKIN_SURFACE_QUADRATIC_SURFACE_H