Signed_mesh_domain_3.h

// Copyright (c) 2009 INRIA Sophia-Antipolis (France).
// 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)     : Stéphane Tayeb (original Labeled_mesh_domain_3.h)
// Author(s)     : Phil Weir (NUMA Engineering modifications)
//
//******************************************************************************
// File Description :
// class Signed_mesh_domain_3. See class description.
//******************************************************************************


#ifndef CGAL_MESH_3_SIGNED_MESH_DOMAIN_3_H
#define CGAL_MESH_3_SIGNED_MESH_DOMAIN_3_H

#if defined(BOOST_MSVC)
#  pragma warning(push)
#  pragma warning(disable:4180) // qualifier applied to function type has no meaning; ignored
#endif

#include <CGAL/Mesh_3/Labeled_mesh_domain_3.h>

namespace CGAL {


template<class Function, class BGT>
class Signed_mesh_domain_3
 : public Mesh_3::Labeled_mesh_domain_3<Function, BGT >
{
public:
  typedef Mesh_3::Labeled_mesh_domain_3<Function, BGT> Base;

  typedef typename Mesh_3::Labeled_mesh_domain_3<Function, BGT>::Subdomain Subdomain;
  typedef typename Mesh_3::Labeled_mesh_domain_3<Function, BGT>::Subdomain_index Subdomain_index;
  typedef typename Mesh_3::Labeled_mesh_domain_3<Function, BGT>::Surface_patch Surface_patch;
  typedef typename Mesh_3::Labeled_mesh_domain_3<Function, BGT>::Segment_3  Segment_3;
  typedef typename Mesh_3::Labeled_mesh_domain_3<Function, BGT>::Ray_3      Ray_3;
  typedef typename Mesh_3::Labeled_mesh_domain_3<Function, BGT>::Line_3     Line_3;
  typedef typename Mesh_3::Labeled_mesh_domain_3<Function, BGT>::Vector_3     Vector_3;
  typedef typename Mesh_3::Labeled_mesh_domain_3<Function, BGT>::Surface_patch_index     Surface_patch_index;
  typedef typename Mesh_3::Labeled_mesh_domain_3<Function, BGT>::Index     Index;
  typedef typename Mesh_3::Labeled_mesh_domain_3<Function, BGT>::Intersection     Intersection;
  typedef typename BGT::Iso_cuboid_3  Iso_cuboid_3;
  typedef typename Base::Sphere_3 Sphere_3;
  typedef typename Base::Point_3 Point_3;
  typedef typename Base::FT FT;

  Signed_mesh_domain_3(const Function& f,
                         const Sphere_3& bounding_sphere,
                         const FT& error_bound = FT(1e-3),
                         CGAL::Random* p_rng = NULL);

  Signed_mesh_domain_3(const Function& f,
                         const Bbox_3& bbox,
                         const FT& error_bound = FT(1e-3),
                         CGAL::Random* p_rng = NULL);

  Signed_mesh_domain_3(const Function& f,
                         const Iso_cuboid_3& bbox,
                         const FT& error_bound = FT(1e-3),
                         CGAL::Random* p_rng = NULL);

  virtual ~Signed_mesh_domain_3()
  {
    if(delete_rng_)
      delete p_rng_;
  }

  struct Construct_initial_points
  {
    Construct_initial_points(const Signed_mesh_domain_3& domain)
      : r_domain_(domain) {}

    template<class OutputIterator>
    OutputIterator operator()(OutputIterator pts, const int n = 12) const;

  private:
    const Signed_mesh_domain_3& r_domain_;
  };

  Construct_initial_points construct_initial_points_object() const
  {
    return Construct_initial_points(*this);
  }

  struct Is_in_domain
  {
    Is_in_domain(const Signed_mesh_domain_3& domain) : r_domain_(domain) {}

    Subdomain operator()(const Point_3& p) const
    {
      // f(p)<=0 means p is outside the domain
      Subdomain_index index = (r_domain_.signed_function_)(p);
      if ( Subdomain_index() >= index )
        return Subdomain();
      else
        return Subdomain(index);
    }
  private:
    const Signed_mesh_domain_3& r_domain_;
  };

  Is_in_domain is_in_domain_object() const { return Is_in_domain(*this); }

  struct Do_intersect_surface
  {
    Do_intersect_surface(const Signed_mesh_domain_3& domain)
      : r_domain_(domain) {}

    Surface_patch operator()(const Segment_3& s) const
    {
      return this->operator()(s.source(), s.target());
    }

    Surface_patch operator()(const Ray_3& r) const
    {
      return clip_to_segment(r);
    }

    Surface_patch operator()(const Line_3& l) const
    {
      return clip_to_segment(l);
    }

  private:
    Surface_patch operator()(const Point_3& a, const Point_3& b) const
    {
      // If f(a) != f(b), then [a,b] intersects some surface. Here we consider
      // [a,b] intersects surface_patch labelled <f(a),f(b)> (or <f(b),f(a)>).
      // It may be false, further rafinement will improve precision
      const Subdomain_index value_a = r_domain_.signed_function_(a);
      const Subdomain_index value_b = r_domain_.signed_function_(b);

      if ( (value_a > 0 || value_b > 0) && value_a != value_b )
        return Surface_patch(r_domain_.make_surface_index(value_a, value_b));
      else
        return Surface_patch();
    }

    template<typename Query>
    Surface_patch clip_to_segment(const Query& query) const
    {
      typename cpp11::result_of<typename BGT::Intersect_3(Query, Iso_cuboid_3)>::type
        clipped = CGAL::intersection(query, r_domain_.bounding_box());

      if(clipped)
#if CGAL_INTERSECTION_VERSION > 1
        if(const Segment_3* s = boost::get<Segment_3>(&*clipped))
          return this->operator()(*s);
#else
        if(const Segment_3* s = object_cast<Segment_3>(&clipped))
          return this->operator()(*s);
#endif
        
      return Surface_patch();
    }

  private:
    const Signed_mesh_domain_3& r_domain_;
  };

  Do_intersect_surface do_intersect_surface_object() const
  {
    return Do_intersect_surface(*this);
  }

  struct Construct_intersection
  {
    Construct_intersection(const Signed_mesh_domain_3& domain)
      : r_domain_(domain) {}

    Intersection operator()(const Segment_3& s) const
    {
#ifndef CGAL_MESH_3_NO_LONGER_CALLS_DO_INTERSECT_3
      CGAL_precondition(r_domain_.do_intersect_surface_object()(s));
#endif // NOT CGAL_MESH_3_NO_LONGER_CALLS_DO_INTERSECT_3
      return this->operator()(s.source(),s.target());
    }

    Intersection operator()(const Ray_3& r) const
    {
      return clip_to_segment(r);
    }

    Intersection operator()(const Line_3& l) const
    {
      return clip_to_segment(l);
    }

  private:
    Intersection operator()(const Point_3& a, const Point_3& b) const
    {
      // Functors
      typename BGT::Compute_squared_distance_3 squared_distance =
                                      BGT().compute_squared_distance_3_object();
      typename BGT::Construct_midpoint_3 midpoint =
                                      BGT().construct_midpoint_3_object();

      // Non const points
      Point_3 p1 = a;
      Point_3 p2 = b;
      Point_3 mid = midpoint(p1, p2);

      // Cannot be const: those values are modified below.
      Subdomain_index value_at_p1 = r_domain_.signed_function_(p1);
      Subdomain_index value_at_p2 = r_domain_.signed_function_(p2);
      Subdomain_index value_at_mid = r_domain_.signed_function_(mid,true);

      // If both extremities are in the same subdomain,
      // there is no intersection.
      // This should not happen...
      if( value_at_p1 == value_at_p2 || (value_at_p1 <= 0 && value_at_p2 <= 0) )
      {
        return Intersection();
      }

      // Construct the surface patch index and index from the values at 'a'
      // and 'b'. Even if the bissection find out a different pair of
      // values, the reported index will be constructed from the initial
      // values.
      const Surface_patch_index sp_index =
        r_domain_.make_surface_index(value_at_p1, value_at_p2);
      const Index index = r_domain_.index_from_surface_patch_index(sp_index);

      // Else lets find a point (by bisection)
      // Bisection ends when the point is near than error bound from surface
      while(true)
      {
        // If the two points are enough close, then we return midpoint
        if ( squared_distance(p1, p2) < r_domain_.squared_error_bound_ )
        {
          CGAL_assertion(value_at_p1 != value_at_p2 && (value_at_p1 > 0 || value_at_p2 > 0));
          return Intersection(mid, index, 2);
        }

        // Else we must go on
        // Here we consider that p1(a) is the source point. Thus, we keep p1 and
        // change p2 if f(p1)!=f(p2).
        // That allows us to find the first intersection from a of [a,b] with
        // a surface.
        if ( value_at_p1 != value_at_mid && (value_at_p1 > 0 || value_at_mid > 0) )
        {
          p2 = mid;
          value_at_p2 = value_at_mid;
        }
        else
        {
          p1 = mid;
          value_at_p1 = value_at_mid;
        }

        mid = midpoint(p1, p2);
        value_at_mid = r_domain_.signed_function_(mid,true);
      }
    }

    template<typename Query>
    Intersection clip_to_segment(const Query& query) const
    {
      typename cpp11::result_of<typename BGT::Intersect_3(Query, Iso_cuboid_3)>::type
        clipped = CGAL::intersection(query, r_domain_.bounding_box());

      if(clipped)
#if CGAL_INTERSECTION_VERSION > 1
        if(const Segment_3* s = boost::get<Segment_3>(&*clipped))
          return this->operator()(*s);
#else
        if(const Segment_3* s = object_cast<Segment_3>(&clipped))
          return this->operator()(*s);
#endif
      
      return Intersection();
    }

  private:
    const Signed_mesh_domain_3& r_domain_;
  };

  Construct_intersection construct_intersection_object() const
  {
    return Construct_intersection(*this);
  }

  Sphere_3 bounding_sphere(const Iso_cuboid_3& bbox) const
  {
    typename BGT::Construct_sphere_3 sphere = BGT().construct_sphere_3_object();
    return sphere((bbox.min)(), (bbox.max)());
  }

protected:
  const Function signed_function_;
  FT squared_error_bound_;

private:
  CGAL::Random* p_rng_;
  bool delete_rng_;

private:
  // Disabled copy constructor & assignment operator
  typedef Signed_mesh_domain_3<Function,BGT> Self;
  Signed_mesh_domain_3(const Self& src);
  Self& operator=(const Self& src);

private:
  Surface_patch_index make_surface_index(const Subdomain_index i,
                                   const Subdomain_index j) const
  {
    if ( i < j ) return Surface_patch_index(i,j);
    else return Surface_patch_index(j,i);
  }

  FT squared_error_bound(const Sphere_3& sphere, const FT& error) const
  {
    typename BGT::Compute_squared_radius_3 squared_radius =
                                    BGT().compute_squared_radius_3_object();
    return squared_radius(sphere)*error*error;
  }

};  // end class Signed_mesh_domain_3


//-------------------------------------------------------
// Method implementation
//-------------------------------------------------------
template<class F, class BGT>
Signed_mesh_domain_3<F,BGT>::Signed_mesh_domain_3(
                       const F& f,
                       const Sphere_3& bounding_sphere,
                       const FT& error_bound,
                       CGAL::Random* p_rng)
: CGAL::Mesh_3::Labeled_mesh_domain_3<F, BGT>(f, bounding_sphere, error_bound)
, signed_function_(f)
, squared_error_bound_(squared_error_bound(bounding_sphere,error_bound))
, p_rng_(p_rng)
, delete_rng_(false)
{
  // TODO : CGAL_ASSERT(0 < f(bounding_sphere.get_center()) ) ?
  if(!p_rng_)
  {
    p_rng_ = new CGAL::Random(0);
    delete_rng_ = true;
  }
}

template<class F, class BGT>
Signed_mesh_domain_3<F,BGT>::Signed_mesh_domain_3(
                       const F& f,
                       const Bbox_3& bbox,
                       const FT& error_bound,
                       CGAL::Random* p_rng)
: CGAL::Mesh_3::Labeled_mesh_domain_3<F, BGT>(f, bbox, error_bound)
, signed_function_(f)
, squared_error_bound_(squared_error_bound(CGAL::Mesh_3::Labeled_mesh_domain_3<F, BGT>::bounding_box(),error_bound))
, p_rng_(p_rng)
, delete_rng_(false)
{
  // TODO : CGAL_ASSERT(0 < f(bounding_sphere.get_center()) ) ?
  if(!p_rng_)
  {
    p_rng_ = new CGAL::Random(0);
    delete_rng_ = true;
  }
}

template<class F, class BGT>
Signed_mesh_domain_3<F,BGT>::Signed_mesh_domain_3(
                       const F& f,
                       const Iso_cuboid_3& bbox,
                       const FT& error_bound,
                       CGAL::Random* p_rng)
: CGAL::Mesh_3::Labeled_mesh_domain_3<F, BGT>(f, bbox, error_bound, p_rng)
, signed_function_(f)
, squared_error_bound_(squared_error_bound(CGAL::Mesh_3::Labeled_mesh_domain_3<F, BGT>::bounding_box(),error_bound))
, p_rng_(p_rng)
, delete_rng_(false)
{
  // TODO : CGAL_ASSERT(0 < f( bbox.get_center()) ) ?
  if(!p_rng_)
  {
    p_rng_ = new CGAL::Random(0);
    delete_rng_ = true;
  }
}

template<class F, class BGT>
template<class OutputIterator>
OutputIterator
Signed_mesh_domain_3<F,BGT>::Construct_initial_points::operator()(
                                                    OutputIterator pts,
                                                    const int nb_points) const
{
  // Create point_iterator on and in bounding_sphere
  typedef Random_points_on_sphere_3<Point_3> Random_points_on_sphere_3;
  typedef Random_points_in_sphere_3<Point_3> Random_points_in_sphere_3;


  const FT squared_radius = BGT().compute_squared_radius_3_object()(
      r_domain_.bounding_sphere(r_domain_.bounding_box()));

  const double radius = std::sqrt(CGAL::to_double(squared_radius));

  CGAL::Random& rng = *(r_domain_.p_rng_);
  Random_points_on_sphere_3 random_point_on_sphere(radius, rng);
  Random_points_in_sphere_3 random_point_in_sphere(radius, rng);

  // Get some functors
  typename BGT::Construct_segment_3 segment_3 =
                              BGT().construct_segment_3_object();
  typename BGT::Construct_vector_3 vector_3 =
                              BGT().construct_vector_3_object();
  typename BGT::Construct_translated_point_3 translate =
                              BGT().construct_translated_point_3_object();
  typename BGT::Construct_center_3 center = BGT().construct_center_3_object();

  // Get translation from origin to sphere center
  Point_3 center_pt = center(r_domain_.bounding_sphere(r_domain_.bounding_box()));
  const Vector_3 sphere_translation = vector_3(CGAL::ORIGIN, center_pt);

  // Create nb_point points
  int n = nb_points;
#ifdef CGAL_MESH_3_VERBOSE
  std::cerr << "construct initial points:\n";
#endif
  while ( 0 != n )
  {
    // Get a random segment
    const Point_3 random_point = translate(*random_point_on_sphere,
                                           sphere_translation);
    const Segment_3 random_seg = segment_3(center_pt, random_point);

    // Add the intersection to the output if it exists
    Surface_patch surface = r_domain_.do_intersect_surface_object()(random_seg);
    if ( surface )
    {
      const Point_3 intersect_pt = CGAL::cpp11::get<0>(
          r_domain_.construct_intersection_object()(random_seg));
      *pts++ = std::make_pair(intersect_pt,
                              r_domain_.index_from_surface_patch_index(*surface));
      --n;

#ifdef CGAL_MESH_3_VERBOSE
      std::cerr << boost::format("\r             \r"
                                 "%1%/%2% initial point(s) found...")
                   % (nb_points - n)
                   % nb_points;
#endif
    }
    else
    {
      // Get a new random point into sphere as center of object
      // It may be necessary if the center of the domain is empty, e.g. torus
      // In general case, it is good for input point dispersion
      ++random_point_in_sphere;
      center_pt = translate(*random_point_in_sphere, sphere_translation);
    }
    ++random_point_on_sphere;
  }

#ifdef CGAL_MESH_3_VERBOSE
  std::cerr << "\n";
#endif
  return pts;
}


}  // end namespace CGAL

#if defined(BOOST_MSVC)
#  pragma warning(pop)
#endif

#endif // CGAL_SIGNED_MESH_DOMAIN_3_H