mirror of https://github.com/CGAL/cgal
Reworked the way periodicity is defined
The input domain does not need to be periodic. It is the domain class that has to handle the periodicity. This is cleaner mathematically and will be more natural for other types of domains. Along the way, the labeled periodic domain is brought up to date with Mesh_3's (bug fixes, null subdomain index, etc.)
This commit is contained in:
Mael Rouxel-Labbé
parent
f222d0fd13
commit
c6f58bce4e
|
|
@ -512,6 +512,7 @@ protected:
|
|||
/// Returns bounding box
|
||||
const Iso_cuboid_3& bounding_box() const { return bbox_; }
|
||||
const Function& labeling_function() const { return function_; }
|
||||
const Null& null_function() const { return null; }
|
||||
FT squared_error_bound_value() const { return squared_error_bound_; }
|
||||
|
||||
private:
|
||||
|
|
|
|||
|
|
@ -51,21 +51,12 @@ using namespace CGAL::parameters;
|
|||
static const int domain_size = 1;
|
||||
static const Iso_cuboid periodic_domain = Iso_cuboid(0, 0, 0, domain_size, domain_size, domain_size);
|
||||
|
||||
// Used to artifically create a periodic function.
|
||||
static const Tr dummy_tr(periodic_domain);
|
||||
Point canonicalize_point(const Point& p)
|
||||
{
|
||||
// Compute the representative of 'p' in the fundamental domain
|
||||
return dummy_tr.canonicalize_point(p);
|
||||
}
|
||||
|
||||
// Implicit function
|
||||
static const FT cx = 0.51, cy = 0.51, cz = 0.5;
|
||||
static const FT scale = 0.9;
|
||||
FT cone_function(const Point& p)
|
||||
{
|
||||
const Point& cp = canonicalize_point(p);
|
||||
const FT x = cp.x(), y = cp.y(), z = cp.z();
|
||||
const FT x = p.x(), y = p.y(), z = p.z();
|
||||
|
||||
if (((x-cx)*(x-cx) + (y-cy)*(y-cy)) > scale * (z-cz)*(z-cz))
|
||||
return 1; // outside
|
||||
|
|
|
|||
|
|
@ -26,6 +26,9 @@
|
|||
|
||||
#include <CGAL/Periodic_3_mesh_3/config.h>
|
||||
|
||||
#include <CGAL/Periodic_3_mesh_triangulation_3.h>
|
||||
#include <CGAL/internal/canonicalize_helper.h>
|
||||
|
||||
#include <CGAL/Labeled_mesh_domain_3.h>
|
||||
#include <CGAL/intersections.h>
|
||||
#include <CGAL/result_of.h>
|
||||
|
|
@ -38,22 +41,31 @@ namespace CGAL
|
|||
/**
|
||||
* \class Labeled_periodic_3_mesh_domain_3
|
||||
*
|
||||
* Function f must take his values into N.
|
||||
* Let p be a Point.
|
||||
* - f(p)=0 means that p is outside domain.
|
||||
* - f(p)=a, a!=0 means that p is inside subdomain a.
|
||||
* Function `f` must be defined over the fundamental domain and take his values into `N`.
|
||||
* Let `p` be a point.
|
||||
* - `f(p)=0` means that p is outside domain.
|
||||
* - `f(p)=a`, `a!=0` means that `p` is inside subdomain `a`.
|
||||
*
|
||||
* Any boundary facet is labelled <a,b>, with a<b, where a and b are the
|
||||
* Any boundary facet is labelled `<a,b>`, with `a<b`, where `a` and `b` are the
|
||||
* tags of its incident subdomain.
|
||||
* Thus, a boundary facet of the domain is labelled <0,b>, where b!=0.
|
||||
* Thus, a boundary facet of the domain is labelled `<0,b>`, where `b!=0`.
|
||||
*/
|
||||
template<class Function, class BGT>
|
||||
// The difference with Mesh_3's is very tiny if 'CGAL_PERIODIC_CANONICALIZE_DUAL_INTERSECTIONS'
|
||||
// is not enabled: it resides in the call to label() which canonicalizes the point
|
||||
// (that is, send them to the fundamental domain) before evaluation.
|
||||
template<class Function, class BGT, class Null_subdomain_index = Default>
|
||||
class Labeled_periodic_3_mesh_domain_3
|
||||
: public Labeled_mesh_domain_3<Function, BGT>
|
||||
: public Labeled_mesh_domain_3<Function, BGT,
|
||||
typename Default::Get<Null_subdomain_index,
|
||||
CGAL::Null_subdomain_index>::type>
|
||||
{
|
||||
public:
|
||||
/// Null subdomain type
|
||||
typedef typename Default::Get<Null_subdomain_index,
|
||||
CGAL::Null_subdomain_index>::type Null;
|
||||
|
||||
/// Base type
|
||||
typedef Labeled_mesh_domain_3<Function, BGT> Base;
|
||||
typedef Labeled_mesh_domain_3<Function, BGT, Null> Base;
|
||||
|
||||
/// Geometric object types
|
||||
typedef typename Base::Point_3 Point_3;
|
||||
|
|
@ -80,13 +92,27 @@ public:
|
|||
typedef typename Base::FT FT;
|
||||
typedef BGT Geom_traits;
|
||||
|
||||
/// Periodic traits used in the canonicalization of the points
|
||||
typedef typename details::Periodic_3_mesh_geom_traits_generator<BGT>::type Periodic_geom_traits;
|
||||
|
||||
Labeled_periodic_3_mesh_domain_3(const Function& f,
|
||||
const Iso_cuboid_3& bbox,
|
||||
const FT& error_bound = FT(1e-6))
|
||||
: Base(f, bbox, error_bound)
|
||||
const FT& error_bound = FT(1e-6),
|
||||
Null null = Null(),
|
||||
CGAL::Random* p_rng = NULL)
|
||||
: Base(f, bbox, error_bound, null, p_rng),
|
||||
pgt(bbox)
|
||||
{ }
|
||||
|
||||
const Iso_cuboid_3& periodic_bounding_box() const { return Base::bounding_box(); }
|
||||
const Periodic_geom_traits& periodic_geom_traits() const { return pgt; }
|
||||
|
||||
Subdomain_index label(const Point_3& p) const {
|
||||
return Base::labeling_function()(P3T3::internal::robust_canonicalize_point(p, pgt));
|
||||
}
|
||||
Subdomain_index label(const Point_3& p, const bool b) const {
|
||||
return Base::labeling_function()(P3T3::internal::robust_canonicalize_point(p, pgt), b);
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns true if the element \ccc{type} intersect properly any of the
|
||||
|
|
@ -127,6 +153,7 @@ public:
|
|||
// [a,b] intersects surface_patch labelled <f(a),f(b)> (or <f(b),f(a)>).
|
||||
// It may be false, further rafinement will improve precision
|
||||
|
||||
// This whole canonicalization of offsets process seems useless... Hiding it behind macros.
|
||||
#ifdef CGAL_PERIODIC_CANONICALIZE_DUAL_INTERSECTIONS
|
||||
Iso_cuboid_3 pbb = r_domain_.periodic_bounding_box();
|
||||
FT dimension [3] = { pbb.xmax()-pbb.xmin(),
|
||||
|
|
@ -146,41 +173,56 @@ public:
|
|||
int o2 [3] = { static_cast<int>(b_t[0] / dimension[0]),
|
||||
static_cast<int>(b_t[1] / dimension[1]),
|
||||
static_cast<int>(b_t[2] / dimension[2]) };
|
||||
#endif
|
||||
|
||||
FT a_min [3] = { a.x(), a.y(), a.z() };
|
||||
FT b_min [3] = { b.x(), b.y(), b.z() };
|
||||
|
||||
#ifdef CGAL_PERIODIC_CANONICALIZE_DUAL_INTERSECTIONS
|
||||
for (unsigned idx = 0; idx < 3; ++idx)
|
||||
{
|
||||
FT offset = dimension[idx] * static_cast<FT>((std::min)(o1[idx], o2[idx]));
|
||||
a_min[idx] -= offset;
|
||||
b_min[idx] -= offset;
|
||||
}
|
||||
|
||||
const Point_3 pa(a_min[0], a_min[1], a_min[2]);
|
||||
const Point_3 pb(b_min[0], b_min[1], b_min[2]);
|
||||
Subdomain_index value_a = r_domain_.label(pa);
|
||||
Subdomain_index value_b = r_domain_.label(pb);
|
||||
#else
|
||||
Subdomain_index value_a = r_domain_.label(a);
|
||||
Subdomain_index value_b = r_domain_.label(b);
|
||||
#endif
|
||||
|
||||
Subdomain_index value_a = r_domain_.labeling_function()(Point_3(a_min[0], a_min[1], a_min[2]));
|
||||
Subdomain_index value_b = r_domain_.labeling_function()(Point_3(b_min[0], b_min[1], b_min[2]));
|
||||
if ( value_a != value_b )
|
||||
return Surface_patch(r_domain_.make_surface_index(value_a, value_b));
|
||||
{
|
||||
if( r_domain_.null_function()(value_a) && r_domain_.null_function()(value_b) )
|
||||
return Surface_patch();
|
||||
else
|
||||
return Surface_patch(r_domain_.make_surface_index(value_a, value_b));
|
||||
}
|
||||
|
||||
#ifdef CGAL_PERIODIC_CANONICALIZE_DUAL_INTERSECTIONS
|
||||
FT a_max [3] = { a.x(), a.y(), a.z() };
|
||||
FT b_max [3] = { b.x(), b.y(), b.z() };
|
||||
|
||||
#ifdef CGAL_PERIODIC_CANONICALIZE_DUAL_INTERSECTIONS
|
||||
for (unsigned idx = 0; idx < 3; ++idx)
|
||||
{
|
||||
FT offset = dimension[idx] * static_cast<FT>((std::max)(o1[idx], o2[idx]));
|
||||
a_max[idx] -= offset;
|
||||
b_max[idx] -= offset;
|
||||
}
|
||||
#endif
|
||||
|
||||
value_a = r_domain_.labeling_function()(Point_3(a_max[0], a_max[1], a_max[2]));
|
||||
value_b = r_domain_.labeling_function()(Point_3(b_max[0], b_max[1], b_max[2]));
|
||||
pa = Point_3(a_max[0], a_max[1], a_max[2]);
|
||||
pb = Point_3(b_max[0], b_max[1], b_max[2]);
|
||||
value_a = r_domain_.label(pa);
|
||||
value_b = r_domain_.label(pb);
|
||||
if ( value_a != value_b )
|
||||
return Surface_patch(r_domain_.make_surface_index(value_a, value_b));
|
||||
{
|
||||
if( r_domain_.null_function()(value_a) && r_domain_.null_function()(value_b) )
|
||||
return Surface_patch();
|
||||
else
|
||||
return Surface_patch(r_domain_.make_surface_index(value_a, value_b));
|
||||
}
|
||||
#endif
|
||||
|
||||
return Surface_patch();
|
||||
}
|
||||
|
|
@ -261,12 +303,10 @@ public:
|
|||
{
|
||||
// Functors
|
||||
typename BGT::Compute_squared_distance_3 squared_distance =
|
||||
BGT().compute_squared_distance_3_object();
|
||||
r_domain_.periodic_geom_traits().compute_squared_distance_3_object();
|
||||
typename BGT::Construct_midpoint_3 midpoint =
|
||||
BGT().construct_midpoint_3_object();
|
||||
r_domain_.periodic_geom_traits().construct_midpoint_3_object();
|
||||
|
||||
// This whole normalization of offsets process seems completely unneeded;
|
||||
// hiding it behind macros.
|
||||
#ifdef CGAL_PERIODIC_CANONICALIZE_DUAL_INTERSECTIONS
|
||||
Iso_cuboid_3 pbb = r_domain_.periodic_bounding_box();
|
||||
FT dimension [3] = { pbb.xmax() - pbb.xmin(),
|
||||
|
|
@ -287,29 +327,30 @@ public:
|
|||
int o2 [3] = { static_cast<int>(b_t[0] / dimension[0]),
|
||||
static_cast<int>(b_t[1] / dimension[1]),
|
||||
static_cast<int>(b_t[2] / dimension[2]) };
|
||||
#endif
|
||||
|
||||
FT a_min [3] = { a.x(), a.y(), a.z() };
|
||||
FT b_min [3] = { b.x(), b.y(), b.z() };
|
||||
|
||||
#ifdef CGAL_PERIODIC_CANONICALIZE_DUAL_INTERSECTIONS
|
||||
for (unsigned idx = 0; idx < 3; ++idx)
|
||||
{
|
||||
FT offset = dimension[idx] * static_cast<FT>((std::min)(o1[idx], o2[idx]));
|
||||
a_min[idx] -= offset;
|
||||
b_min[idx] -= offset;
|
||||
}
|
||||
#endif
|
||||
|
||||
// Non const points
|
||||
Point_3 p1(a_min[0], a_min[1], a_min[2]);
|
||||
Point_3 p2(b_min[0], b_min[1], b_min[2]);
|
||||
#else
|
||||
Point_3 p1 = a;
|
||||
Point_3 p2 = b;
|
||||
#endif
|
||||
|
||||
Point_3 mid = midpoint(p1, p2);
|
||||
|
||||
// Cannot be const: those values are modified below.
|
||||
Subdomain_index value_at_p1 = r_domain_.labeling_function()(p1);
|
||||
Subdomain_index value_at_p2 = r_domain_.labeling_function()(p2);
|
||||
Subdomain_index value_at_mid = r_domain_.labeling_function()(mid, true);
|
||||
Subdomain_index value_at_p1 = r_domain_.label(p1);
|
||||
Subdomain_index value_at_p2 = r_domain_.label(p2);
|
||||
Subdomain_index value_at_mid = r_domain_.label(mid, true);
|
||||
|
||||
// If both extremities are in the same subdomain, then there is no intersection.
|
||||
// This should not happen...
|
||||
|
|
@ -324,26 +365,23 @@ public:
|
|||
a_max[idx] -= offset;
|
||||
b_max[idx] -= offset;
|
||||
}
|
||||
|
||||
p1 = Point_3(a_max[0], a_max[1], a_max[2]);
|
||||
p2 = Point_3(b_max[0], b_max[1], b_max[2]);
|
||||
mid = midpoint(p1, p2);
|
||||
#endif
|
||||
|
||||
value_at_p1 = r_domain_.labeling_function()(p1);
|
||||
value_at_p2 = r_domain_.labeling_function()(p2);
|
||||
value_at_mid = r_domain_.labeling_function()(mid, true);
|
||||
value_at_p1 = r_domain_.label(p1);
|
||||
value_at_p2 = r_domain_.label(p2);
|
||||
value_at_mid = r_domain_.label(mid, true);
|
||||
|
||||
if( value_at_p1 == value_at_p2 )
|
||||
#endif
|
||||
return Intersection();
|
||||
}
|
||||
|
||||
// Construct the surface patch index and index from the values at 'a'
|
||||
// and 'b'. Even if the bissection finds 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);
|
||||
if( r_domain_.null_function()(value_at_p1) && r_domain_.null_function()(value_at_p2) ) {
|
||||
return Intersection();
|
||||
}
|
||||
|
||||
// Else lets find a point (by bisection)
|
||||
// Bisection ends when the point is closer to the surface than the error bound
|
||||
|
|
@ -352,7 +390,13 @@ public:
|
|||
// If the two points are sufficiently close, then we return the midpoint
|
||||
if ( squared_distance(p1, p2) < r_domain_.squared_error_bound_value() )
|
||||
{
|
||||
CGAL_assertion(value_at_p1 != value_at_p2);
|
||||
CGAL_assertion(value_at_p1 != value_at_p2 &&
|
||||
! ( r_domain_.null_function()(value_at_p1) && r_domain_.null_function()(value_at_p2) ) );
|
||||
|
||||
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);
|
||||
|
||||
return Intersection(mid, index, 2);
|
||||
}
|
||||
|
||||
|
|
@ -361,7 +405,8 @@ public:
|
|||
// 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 )
|
||||
if ( value_at_p1 != value_at_mid &&
|
||||
! ( r_domain_.null_function()(value_at_p1) && r_domain_.null_function()(value_at_mid) ) )
|
||||
{
|
||||
p2 = mid;
|
||||
value_at_p2 = value_at_mid;
|
||||
|
|
@ -373,7 +418,7 @@ public:
|
|||
}
|
||||
|
||||
mid = midpoint(p1, p2);
|
||||
value_at_mid = r_domain_.labeling_function()(mid, true);
|
||||
value_at_mid = r_domain_.label(mid, true);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -405,6 +450,14 @@ public:
|
|||
return Construct_intersection(*this);
|
||||
}
|
||||
|
||||
private:
|
||||
// The domain must know the periodic fundamental domain, which is in the traits.
|
||||
//
|
||||
// Note that the fundamental domain is also known through 'bounding_box()', but
|
||||
// a full geometric traits class is needed for 'construct_point()'
|
||||
// in canonicalization functions.
|
||||
Periodic_geom_traits pgt;
|
||||
|
||||
private:
|
||||
// Disabled copy constructor & assignment operator
|
||||
Labeled_periodic_3_mesh_domain_3(const Labeled_periodic_3_mesh_domain_3&);
|
||||
|
|
|
|||
|
|
@ -32,6 +32,7 @@
|
|||
#include <CGAL/Kernel_traits.h>
|
||||
#include <CGAL/Robust_weighted_circumcenter_filtered_traits_3.h>
|
||||
#include <CGAL/internal/Robust_periodic_weighted_circumcenter_traits_3.h>
|
||||
#include <CGAL/internal/canonicalize_helper.h>
|
||||
|
||||
// periodic triangulations
|
||||
#include <CGAL/Periodic_3_regular_triangulation_traits_3.h>
|
||||
|
|
@ -194,152 +195,19 @@ public:
|
|||
|
||||
void set_domain(const Iso_cuboid& domain)
|
||||
{
|
||||
const FT eps = std::numeric_limits<FT>::epsilon(); // @tmp
|
||||
std::cout << "epsilon: " << eps << std::endl;
|
||||
Base::set_domain(domain);
|
||||
}
|
||||
|
||||
bool is_point_too_close_to_border(const Periodic_bare_point& pbp) const
|
||||
{
|
||||
Bare_point p = construct_point(pbp);
|
||||
const FT px = p.x();
|
||||
const FT py = p.y();
|
||||
const FT pz = p.z();
|
||||
const FT dxm = domain().xmin();
|
||||
const FT dym = domain().ymin();
|
||||
const FT dzm = domain().zmin();
|
||||
const FT dxM = domain().xmax();
|
||||
const FT dyM = domain().ymax();
|
||||
const FT dzM = domain().zmax();
|
||||
|
||||
// simply comparing to FT::epsilon() is probably not completely satisfactory
|
||||
const FT eps = std::numeric_limits<FT>::epsilon();
|
||||
|
||||
FT diff = CGAL::abs(px - dxm);
|
||||
if(diff < eps && diff > 0) return true;
|
||||
diff = CGAL::abs(px - dxM);
|
||||
if(diff < eps && diff > 0) return true;
|
||||
diff = CGAL::abs(py - dym);
|
||||
if(diff < eps && diff > 0) return true;
|
||||
diff = CGAL::abs(py - dyM);
|
||||
if(diff < eps && diff > 0) return true;
|
||||
diff = CGAL::abs(pz - dzm);
|
||||
if(diff < eps && diff > 0) return true;
|
||||
diff = CGAL::abs(pz - dzM);
|
||||
if(diff < eps && diff > 0) return true;
|
||||
return false;
|
||||
}
|
||||
|
||||
Bare_point snap_to_domain_border(const Bare_point& p) const
|
||||
{
|
||||
const FT px = p.x();
|
||||
const FT py = p.y();
|
||||
const FT pz = p.z();
|
||||
const FT dxm = domain().xmin();
|
||||
const FT dym = domain().ymin();
|
||||
const FT dzm = domain().zmin();
|
||||
const FT dxM = domain().xmax();
|
||||
const FT dyM = domain().ymax();
|
||||
const FT dzM = domain().zmax();
|
||||
FT sx = px, sy = py, sz = pz;
|
||||
|
||||
// simply comparing to FT::epsilon() is probably not completely satisfactory
|
||||
const FT eps = std::numeric_limits<FT>::epsilon();
|
||||
|
||||
if(CGAL::abs(px - dxm) < eps) sx = domain().xmin();
|
||||
if(CGAL::abs(px - dxM) < eps) sx = domain().xmax();
|
||||
if(CGAL::abs(py - dym) < eps) sy = domain().ymin();
|
||||
if(CGAL::abs(py - dyM) < eps) sy = domain().ymax();
|
||||
if(CGAL::abs(pz - dzm) < eps) sz = domain().zmin();
|
||||
if(CGAL::abs(pz - dzM) < eps) sz = domain().zmax();
|
||||
|
||||
std::cout << "epsilon: " << eps << std::endl;
|
||||
std::cout.precision(20);
|
||||
std::cout << "snapped " << p << " to " << sx << " " << sy << " " << sz << std::endl;
|
||||
|
||||
return geom_traits().construct_point_3_object()(sx, sy, sz);
|
||||
}
|
||||
|
||||
Weighted_point snap_to_domain_border(const Weighted_point& p) const
|
||||
{
|
||||
typename Geom_traits::Compute_weight_3 cw = geom_traits().compute_weight_3_object();
|
||||
|
||||
const Bare_point snapped_p =
|
||||
snap_to_domain_border(geom_traits().construct_point_3_object()(p));
|
||||
|
||||
return geom_traits().construct_weighted_point_3_object()(snapped_p, cw(p));
|
||||
}
|
||||
|
||||
/// transform a bare point (living anywhere in space) into the canonical
|
||||
/// instance of the same bare point that lives inside the base domain
|
||||
Bare_point robust_canonicalize_point(const Bare_point& p) const
|
||||
{
|
||||
bool should_snap = false;
|
||||
Periodic_bare_point pbp = construct_periodic_point(p, should_snap);
|
||||
|
||||
if(!should_snap)
|
||||
{
|
||||
// Even if there is no issue while constructing the canonical point,
|
||||
// snap the point if it's too close to a border of the domain
|
||||
should_snap = is_point_too_close_to_border(pbp);
|
||||
}
|
||||
|
||||
if(should_snap)
|
||||
{
|
||||
Bare_point sp = snap_to_domain_border(p);
|
||||
|
||||
// might have snapped to a 'max' of the domain, which is not in the domain
|
||||
// note: we could snap to 'min' all the time in 'snap_to_domain_border'
|
||||
// but this is clearer like that (and costs very little since we should
|
||||
// not have to use exact computations too often)
|
||||
return canonicalize_point(sp);
|
||||
}
|
||||
|
||||
Bare_point canonical_p = construct_point(pbp);
|
||||
CGAL_postcondition( !(canonical_p.x() < domain().xmin()) &&
|
||||
(canonical_p.x() < domain().xmax()) );
|
||||
CGAL_postcondition( !(canonical_p.y() < domain().ymin()) &&
|
||||
(canonical_p.y() < domain().ymax()) );
|
||||
CGAL_postcondition( !(canonical_p.z() < domain().zmin()) &&
|
||||
(canonical_p.z() < domain().zmax()) );
|
||||
|
||||
return canonical_p;
|
||||
}
|
||||
|
||||
Bare_point canonicalize_point(const Bare_point& p) const
|
||||
{
|
||||
if(p.x() >= domain().xmin() && p.x() < domain().xmax() &&
|
||||
p.y() >= domain().ymin() && p.y() < domain().ymax() &&
|
||||
p.z() >= domain().zmin() && p.z() < domain().zmax())
|
||||
return p;
|
||||
|
||||
return robust_canonicalize_point(p);
|
||||
}
|
||||
|
||||
/// transform a weighted point (living anywhere in space) into the canonical
|
||||
/// instance of the same weighted point that lives inside the base domain
|
||||
Weighted_point robust_canonicalize_point(const Weighted_point& p) const
|
||||
{
|
||||
typename Geom_traits::Compute_weight_3 cw = geom_traits().compute_weight_3_object();
|
||||
typename Geom_traits::Construct_point_3 cp = geom_traits().construct_point_3_object();
|
||||
typename Geom_traits::Construct_weighted_point_3 cwp = geom_traits().construct_weighted_point_3_object();
|
||||
|
||||
const Bare_point& bp = cp(p);
|
||||
Bare_point canonical_point = robust_canonicalize_point(bp);
|
||||
|
||||
return cwp(canonical_point, cw(p));
|
||||
return P3T3::internal::robust_canonicalize_point(p, geom_traits());
|
||||
}
|
||||
|
||||
// @todo it might be dangerous to call robust_canonicalize without also changing
|
||||
// <p, offset> = construct_periodic_point(p) (lack of consistency in the result)
|
||||
Weighted_point canonicalize_point(const Weighted_point& p) const
|
||||
{
|
||||
if(p.x() >= domain().xmin() && p.x() < domain().xmax() &&
|
||||
p.y() >= domain().ymin() && p.y() < domain().ymax() &&
|
||||
p.z() >= domain().zmin() && p.z() < domain().zmax())
|
||||
return p;
|
||||
|
||||
// @todo it might be dangerous to call robust_canonicalize without also changing
|
||||
// <p, offset> = construct_periodic_point(p) (lack of consistency in the result)
|
||||
return robust_canonicalize_point(p);
|
||||
return P3T3::internal::robust_canonicalize_point(p, geom_traits());
|
||||
}
|
||||
|
||||
Triangle triangle(const Facet& f) const
|
||||
|
|
|
|||
|
|
@ -37,16 +37,16 @@ typedef CGAL::Mesh_criteria_3<Tr> Periodic_mesh_criteria;
|
|||
// To avoid verbose function and named parameters call
|
||||
using namespace CGAL::parameters;
|
||||
|
||||
static Tr dummy_tr; // default constructs a iso_cuboid (0,0,0,1,1,1)
|
||||
|
||||
// Function
|
||||
FT sphere(const Point& p)
|
||||
{
|
||||
assert(p.x() >= 0 && p.y() >= 0 && p.z() >= 0 && p.x() < 1 && p.y() < 1 && p.z() < 1);
|
||||
return CGAL::squared_distance(p, Point(0.5, 0.5, 0.5))-0.2;
|
||||
}
|
||||
|
||||
FT schwarz_p(const Point& p)
|
||||
{
|
||||
assert(p.x() >= 0 && p.y() >= 0 && p.z() >= 0 && p.x() < 1 && p.y() < 1 && p.z() < 1);
|
||||
const FT x2 = std::cos( p.x() * 2 * CGAL_PI ),
|
||||
y2 = std::cos( p.y() * 2 * CGAL_PI ),
|
||||
z2 = std::cos( p.z() * 2 * CGAL_PI );
|
||||
|
|
@ -56,11 +56,13 @@ FT schwarz_p(const Point& p)
|
|||
// not triply periodic
|
||||
FT scherk(const Point& p)
|
||||
{
|
||||
assert(p.x() >= 0 && p.y() >= 0 && p.z() >= 0 && p.x() < 1 && p.y() < 1 && p.z() < 1);
|
||||
return std::exp( p.z() ) * std::cos( p.y() * 2 * CGAL_PI) - std::cos( p.x() * 2 * CGAL_PI );
|
||||
}
|
||||
|
||||
// Triply Implicit Periodic Functions for meshing
|
||||
FT schwarz_p_transl (const Point& p) {
|
||||
assert(p.x() >= 0 && p.y() >= 0 && p.z() >= 0 && p.x() < 1 && p.y() < 1 && p.z() < 1);
|
||||
const FT x2 = std::cos(p.x() * 2 * CGAL_PI + CGAL_PI / 2.0),
|
||||
y2 = std::cos(p.y() * 2 * CGAL_PI + CGAL_PI / 2.0),
|
||||
z2 = std::cos(p.z() * 2 * CGAL_PI + CGAL_PI / 2.0);
|
||||
|
|
@ -68,6 +70,7 @@ FT schwarz_p_transl (const Point& p) {
|
|||
}
|
||||
|
||||
FT gyroid (const Point& p) {
|
||||
assert(p.x() >= 0 && p.y() >= 0 && p.z() >= 0 && p.x() < 1 && p.y() < 1 && p.z() < 1);
|
||||
const FT cx = std::cos(p.x() * 2 * CGAL_PI),
|
||||
cy = std::cos(p.y() * 2 * CGAL_PI),
|
||||
cz = std::cos(p.z() * 2 * CGAL_PI);
|
||||
|
|
@ -78,6 +81,7 @@ FT gyroid (const Point& p) {
|
|||
}
|
||||
|
||||
FT diamond (const Point& p) {
|
||||
assert(p.x() >= 0 && p.y() >= 0 && p.z() >= 0 && p.x() < 1 && p.y() < 1 && p.z() < 1);
|
||||
const FT cx = std::cos(p.x() * 2 * CGAL_PI),
|
||||
cy = std::cos(p.y() * 2 * CGAL_PI),
|
||||
cz = std::cos(p.z() * 2 * CGAL_PI);
|
||||
|
|
@ -88,6 +92,7 @@ FT diamond (const Point& p) {
|
|||
}
|
||||
|
||||
FT double_p (const Point& p) {
|
||||
assert(p.x() >= 0 && p.y() >= 0 && p.z() >= 0 && p.x() < 1 && p.y() < 1 && p.z() < 1);
|
||||
const FT cx = std::cos(p.x() * 2 * CGAL_PI),
|
||||
cy = std::cos(p.y() * 2 * CGAL_PI),
|
||||
cz = std::cos(p.z() * 2 * CGAL_PI);
|
||||
|
|
@ -98,6 +103,7 @@ FT double_p (const Point& p) {
|
|||
}
|
||||
|
||||
FT G_prime (const Point& p) {
|
||||
assert(p.x() >= 0 && p.y() >= 0 && p.z() >= 0 && p.x() < 1 && p.y() < 1 && p.z() < 1);
|
||||
const FT cx = std::cos(p.x() * 2 * CGAL_PI),
|
||||
cy = std::cos(p.y() * 2 * CGAL_PI),
|
||||
cz = std::cos(p.z() * 2 * CGAL_PI);
|
||||
|
|
@ -115,6 +121,7 @@ FT G_prime (const Point& p) {
|
|||
}
|
||||
|
||||
FT lidinoid (const Point& p) {
|
||||
assert(p.x() >= 0 && p.y() >= 0 && p.z() >= 0 && p.x() < 1 && p.y() < 1 && p.z() < 1);
|
||||
const FT cx = std::cos(p.x() * 2 * CGAL_PI),
|
||||
cy = std::cos(p.y() * 2 * CGAL_PI),
|
||||
cz = std::cos(p.z() * 2 * CGAL_PI);
|
||||
|
|
@ -132,6 +139,7 @@ FT lidinoid (const Point& p) {
|
|||
}
|
||||
|
||||
FT D_prime (const Point& p) {
|
||||
assert(p.x() >= 0 && p.y() >= 0 && p.z() >= 0 && p.x() < 1 && p.y() < 1 && p.z() < 1);
|
||||
const FT cx = std::cos(p.x() * 2 * CGAL_PI),
|
||||
cy = std::cos(p.y() * 2 * CGAL_PI),
|
||||
cz = std::cos(p.z() * 2 * CGAL_PI);
|
||||
|
|
@ -146,6 +154,7 @@ FT D_prime (const Point& p) {
|
|||
}
|
||||
|
||||
FT split_p (const Point& p) {
|
||||
assert(p.x() >= 0 && p.y() >= 0 && p.z() >= 0 && p.x() < 1 && p.y() < 1 && p.z() < 1);
|
||||
const FT cx = std::cos(p.x() * 2 * CGAL_PI),
|
||||
cy = std::cos(p.y() * 2 * CGAL_PI),
|
||||
cz = std::cos(p.z() * 2 * CGAL_PI);
|
||||
|
|
@ -163,22 +172,11 @@ FT split_p (const Point& p) {
|
|||
- 0.4 * (cx + cy + cz);
|
||||
}
|
||||
|
||||
// simple implementation
|
||||
Point canonicalize_point(const Point& p)
|
||||
{
|
||||
Point canonical_p = dummy_tr.canonicalize_point(p);
|
||||
|
||||
CGAL_postcondition( !(canonical_p.x() < 0) && (canonical_p.x() < 1) );
|
||||
CGAL_postcondition( !(canonical_p.y() < 0) && (canonical_p.y() < 1) );
|
||||
CGAL_postcondition( !(canonical_p.z() < 0) && (canonical_p.z() < 1) );
|
||||
return canonical_p;
|
||||
}
|
||||
|
||||
FT cylinder(const Point& p)
|
||||
{
|
||||
Point p_ = canonicalize_point(p);
|
||||
const FT x2 = (p_.y() - 0.5) * (p_.y() - 0.5),
|
||||
y2 = (p_.z() - 0.5) * (p_.z() - 0.5);
|
||||
assert(p.x() >= 0 && p.y() >= 0 && p.z() >= 0 && p.x() < 1 && p.y() < 1 && p.z() < 1);
|
||||
const FT x2 = (p.y() - 0.5) * (p.y() - 0.5),
|
||||
y2 = (p.z() - 0.5) * (p.z() - 0.5);
|
||||
return x2 + y2 - 0.1;
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -45,25 +45,11 @@ typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteri
|
|||
// To avoid verbose function and named parameters call
|
||||
using namespace CGAL::parameters;
|
||||
|
||||
// Function
|
||||
static Tr dummy_tr; // default constructs a iso_cuboid (0,0,0,1,1,1)
|
||||
|
||||
Point canonicalize_point(const Point& p)
|
||||
{
|
||||
Point canonical_p = dummy_tr.robust_canonicalize_point(p);
|
||||
|
||||
CGAL_postcondition( !(canonical_p.x() < 0) && (canonical_p.x() < 1) );
|
||||
CGAL_postcondition( !(canonical_p.y() < 0) && (canonical_p.y() < 1) );
|
||||
CGAL_postcondition( !(canonical_p.z() < 0) && (canonical_p.z() < 1) );
|
||||
return canonical_p;
|
||||
}
|
||||
|
||||
////////////////////////////////////////////////////////////////////////////////
|
||||
static const FT r = 0.6;
|
||||
|
||||
FT sphere_function(const Point& p)
|
||||
{
|
||||
return CGAL::squared_distance(canonicalize_point(p), Point(0.5, 0.5, 0.5)) - r*r;
|
||||
return CGAL::squared_distance(p, Point(0.5, 0.5, 0.5)) - r*r;
|
||||
}
|
||||
|
||||
void fill_sphere_polylines(Polylines& polylines)
|
||||
|
|
|
|||
|
|
@ -29,7 +29,6 @@
|
|||
#include <CGAL/license/Periodic_3_triangulation_3.h>
|
||||
|
||||
#include <CGAL/basic.h>
|
||||
#include <CGAL/Exact_kernel_selector.h>
|
||||
|
||||
#include <boost/tuple/tuple.hpp>
|
||||
#include <boost/random/linear_congruential.hpp>
|
||||
|
|
@ -44,8 +43,10 @@
|
|||
#include <CGAL/Triangulation_cell_base_3.h>
|
||||
#include <CGAL/Triangulation_vertex_base_3.h>
|
||||
#include <CGAL/triangulation_assertions.h>
|
||||
#include <CGAL/internal/canonicalize_helper.h>
|
||||
|
||||
#include <CGAL/array.h>
|
||||
#include <CGAL/internal/Exact_type_selector.h>
|
||||
#include <CGAL/NT_converter.h>
|
||||
#include <CGAL/Unique_hash_map.h>
|
||||
#include <CGAL/use.h>
|
||||
|
|
@ -614,159 +615,13 @@ public:
|
|||
return construct_point(pp.first, pp.second);
|
||||
}
|
||||
|
||||
Periodic_point_3 construct_periodic_point_exact(const Point_3& p) const
|
||||
{
|
||||
typedef typename Geometric_traits::Kernel K;
|
||||
typedef typename Exact_kernel_selector<K>::Exact_kernel EK;
|
||||
typedef typename Exact_kernel_selector<K>::C2E C2E;
|
||||
|
||||
C2E to_exact;
|
||||
|
||||
typedef Periodic_3_triangulation_traits_3<EK> Exact_traits;
|
||||
Exact_traits etraits(to_exact(domain()));
|
||||
|
||||
Offset transl(0, 0, 0);
|
||||
typename EK::Point_3 ep = to_exact(p);
|
||||
typename EK::Point_3 dp;
|
||||
|
||||
const typename EK::Iso_cuboid_3& dom = etraits.get_domain();
|
||||
|
||||
while(true) /* while not in */
|
||||
{
|
||||
dp = etraits.construct_point_3_object()(ep, transl);
|
||||
|
||||
if(dp.x() < dom.xmin())
|
||||
transl.x() += 1;
|
||||
else if(dp.y() < dom.ymin())
|
||||
transl.y() += 1;
|
||||
else if(dp.z() < dom.zmin())
|
||||
transl.z() += 1;
|
||||
else if(!(dp.x() < dom.xmax()))
|
||||
transl.x() -= 1;
|
||||
else if(!(dp.y() < dom.ymax()))
|
||||
transl.y() -= 1;
|
||||
else if(!(dp.z() < dom.zmax()))
|
||||
transl.z() -= 1;
|
||||
else
|
||||
break;
|
||||
}
|
||||
|
||||
Periodic_point_3 pp(std::make_pair(p, transl));
|
||||
|
||||
return pp;
|
||||
}
|
||||
|
||||
// Given a point `p` in space, compute its offset `o` with respect
|
||||
// to the canonical domain and returns (p, o)
|
||||
Periodic_point_3 construct_periodic_point(const Point_3& p,
|
||||
bool& encountered_issue) const
|
||||
{
|
||||
// Check if p lies within the domain. If not, translate.
|
||||
const Iso_cuboid& dom = domain();
|
||||
if(!(p.x() < dom.xmin()) && p.x() < dom.xmax() &&
|
||||
!(p.y() < dom.ymin()) && p.y() < dom.ymax() &&
|
||||
!(p.z() < dom.zmin()) && p.z() < dom.zmax())
|
||||
return std::make_pair(p, Offset());
|
||||
|
||||
// Numerical approximations might create inconsistencies between the constructions
|
||||
// and the comparisons. For example in a cubic domain of size 2:
|
||||
// 1. initial point: P(2+1e-17, 0, 0)
|
||||
// 2. the function computes an offset(1, 0, 0),
|
||||
// 3. P + (-1, 0, 0) * domain_size constructs Q(-1e-17, 0, 0) // numerical approximation
|
||||
// 4. the function computes an offset of (-1, 0, 0)
|
||||
// 5. Q + (1, 0, 0) * domain_size constructs (2+1e-17, 0, 0) (that is P)
|
||||
// And the function is looping...
|
||||
//
|
||||
// If this is happening the 'Last_change' enum will break this infinite
|
||||
// loop and return the wrong point and the 'encountered_issue' bool will be
|
||||
// set to 'true'. An exact version of this function should then be called.
|
||||
|
||||
enum Last_change {
|
||||
NO_LAST_CHANGE,
|
||||
INCREASED_X, DECREASED_X, INCREASED_Y, DECREASED_Y, INCREASED_Z, DECREASED_Z
|
||||
};
|
||||
|
||||
Last_change lc = NO_LAST_CHANGE;
|
||||
bool in = false;
|
||||
|
||||
Offset transl(0, 0, 0);
|
||||
Point_3 dp;
|
||||
|
||||
while(!in)
|
||||
{
|
||||
dp = construct_point(std::make_pair(p, transl));
|
||||
|
||||
if(dp.x() < dom.xmin())
|
||||
{
|
||||
if(lc == DECREASED_X) // stuck in a loop
|
||||
break;
|
||||
|
||||
lc = INCREASED_X;
|
||||
transl.x() += 1;
|
||||
}
|
||||
else if(dp.y() < dom.ymin())
|
||||
{
|
||||
if(lc == DECREASED_Y) // stuck in a loop
|
||||
break;
|
||||
|
||||
lc = INCREASED_Y;
|
||||
transl.y() += 1;
|
||||
}
|
||||
else if(dp.z() < dom.zmin())
|
||||
{
|
||||
if(lc == DECREASED_Z) // stuck in a loop
|
||||
break;
|
||||
|
||||
lc = INCREASED_Z;
|
||||
transl.z() += 1;
|
||||
}
|
||||
else if(!(dp.x() < dom.xmax()))
|
||||
{
|
||||
if(lc == INCREASED_X) // stuck in a loop
|
||||
break;
|
||||
|
||||
lc = DECREASED_X;
|
||||
transl.x() -= 1;
|
||||
}
|
||||
else if(!(dp.y() < dom.ymax()))
|
||||
{
|
||||
if(lc == INCREASED_Y) // stuck in a loop
|
||||
break;
|
||||
|
||||
lc = DECREASED_Y;
|
||||
transl.y() -= 1;
|
||||
}
|
||||
else if(!(dp.z() < dom.zmax()))
|
||||
{
|
||||
if(lc == INCREASED_Z) // stuck in a loop
|
||||
break;
|
||||
|
||||
lc = DECREASED_Z;
|
||||
transl.z() -= 1;
|
||||
}
|
||||
else
|
||||
{
|
||||
in = true;
|
||||
}
|
||||
}
|
||||
|
||||
Periodic_point_3 pp(std::make_pair(p, transl));
|
||||
|
||||
if(dp.x() < dom.xmin() || !(dp.x() < dom.xmax()) ||
|
||||
dp.y() < dom.ymin() || !(dp.y() < dom.ymax()) ||
|
||||
dp.z() < dom.zmin() || !(dp.z() < dom.zmax()))
|
||||
{
|
||||
encountered_issue = true;
|
||||
pp = construct_periodic_point_exact(p);
|
||||
}
|
||||
|
||||
return pp;
|
||||
}
|
||||
|
||||
Periodic_point_3 construct_periodic_point(const Point_3& p) const
|
||||
{
|
||||
bool useless = false;
|
||||
return construct_periodic_point(p, useless);
|
||||
|
||||
// The function is a different file to be able to be used where there is
|
||||
// no triangulation (namely, the domains of Periodic_3_mesh_3).
|
||||
return ::CGAL::P3T3::internal::construct_periodic_point(p, useless, geom_traits());
|
||||
}
|
||||
|
||||
// ---------------------------------------------------------------------------
|
||||
|
|
|
|||
|
|
@ -0,0 +1,396 @@
|
|||
// Copyright (c) 2017 GeometryFactory (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$
|
||||
// SPDX-License-Identifier: GPL-3.0+
|
||||
//
|
||||
// Author(s) : Mael Rouxel-Labbé
|
||||
|
||||
// The purpose of this file is to implement a function 'canonicalize_point()'
|
||||
// which transforms a point into its canonical representative (the representative
|
||||
// that belongs to the domain). To improve robustness, the function is allowed
|
||||
// to modify slightly the position of the point, snapping it to the domain if
|
||||
// it is epsilon-close (and outside).
|
||||
//
|
||||
// Although 'canoncalize_point()' is used by Periodic_3_mesh_3, this file is here
|
||||
// (in Periodic_3_triangulation_3) because of 'construct_periodic_point()',
|
||||
// which is a function used in P3T3.h and also needed by 'canonicalize_point()'.
|
||||
// However, P3M3 needs 'canoncalize_point()' without having access to a triangulation
|
||||
// and to avoid duplicating it, the function is here.
|
||||
|
||||
// Geom_traits must be a model of the concept 'P3T3Traits' for 'construct_periodic_point()'.
|
||||
// Geom_traits must be a model of the concept 'P3T3RegularTraits' for everything else.
|
||||
|
||||
#ifndef CGAL_PERIODIC_3_TRIANGULATION_3_CANONICALIZE_HELPER_H
|
||||
#define CGAL_PERIODIC_3_TRIANGULATION_3_CANONICALIZE_HELPER_H
|
||||
|
||||
#include <CGAL/license/Periodic_3_triangulation_3.h>
|
||||
|
||||
#include <CGAL/Periodic_3_triangulation_traits_3.h>
|
||||
|
||||
#include <CGAL/Exact_kernel_selector.h>
|
||||
|
||||
#include <utility>
|
||||
|
||||
namespace CGAL {
|
||||
namespace P3T3 { // can't name it Periodic_3_triangulation_3 because it's already a class...
|
||||
namespace internal {
|
||||
|
||||
template <typename Gt_>
|
||||
std::pair<typename Gt_::Point_3, typename Gt_::Periodic_3_offset_3>
|
||||
construct_periodic_point_exact(const typename Gt_::Point_3& p,
|
||||
const Gt_& gt)
|
||||
{
|
||||
typedef Gt_ Geom_traits;
|
||||
typedef typename Geom_traits::Periodic_3_offset_3 Offset;
|
||||
typedef typename Geom_traits::Iso_cuboid_3 Iso_cuboid;
|
||||
|
||||
const Iso_cuboid& domain = gt.get_domain();
|
||||
|
||||
typedef typename Geom_traits::Kernel K;
|
||||
typedef typename Exact_kernel_selector<K>::Exact_kernel EK;
|
||||
typedef typename Exact_kernel_selector<K>::C2E C2E;
|
||||
|
||||
C2E to_exact;
|
||||
|
||||
typedef Periodic_3_triangulation_traits_3<EK> Exact_traits;
|
||||
Exact_traits etraits(to_exact(domain));
|
||||
|
||||
Offset transl(0, 0, 0);
|
||||
typename EK::Point_3 ep = to_exact(p);
|
||||
typename EK::Point_3 dp;
|
||||
|
||||
const typename EK::Iso_cuboid_3& exact_domain = etraits.get_domain();
|
||||
|
||||
while(true) /* while not in */
|
||||
{
|
||||
dp = etraits.construct_point_3_object()(ep, transl);
|
||||
|
||||
if(dp.x() < exact_domain.xmin())
|
||||
transl.x() += 1;
|
||||
else if(dp.y() < exact_domain.ymin())
|
||||
transl.y() += 1;
|
||||
else if(dp.z() < exact_domain.zmin())
|
||||
transl.z() += 1;
|
||||
else if(!(dp.x() < exact_domain.xmax()))
|
||||
transl.x() -= 1;
|
||||
else if(!(dp.y() < exact_domain.ymax()))
|
||||
transl.y() -= 1;
|
||||
else if(!(dp.z() < exact_domain.zmax()))
|
||||
transl.z() -= 1;
|
||||
else
|
||||
break;
|
||||
}
|
||||
|
||||
return std::make_pair(p, transl);
|
||||
}
|
||||
|
||||
// Given a point `p` in space, compute its offset `o` with respect
|
||||
// to the canonical domain and returns (p, o)
|
||||
template <typename Gt_>
|
||||
std::pair<typename Gt_::Point_3, typename Gt_::Periodic_3_offset_3>
|
||||
construct_periodic_point(const typename Gt_::Point_3& p, bool& encountered_issue, const Gt_& gt)
|
||||
{
|
||||
typedef Gt_ Geom_traits;
|
||||
typedef typename Geom_traits::Point_3 Point;
|
||||
typedef typename Geom_traits::Periodic_3_offset_3 Offset;
|
||||
typedef typename Geom_traits::Iso_cuboid_3 Iso_cuboid;
|
||||
|
||||
const Iso_cuboid& domain = gt.get_domain();
|
||||
|
||||
// Check if p lies within the domain. If not, translate.
|
||||
if(!(p.x() < domain.xmin()) && p.x() < domain.xmax() &&
|
||||
!(p.y() < domain.ymin()) && p.y() < domain.ymax() &&
|
||||
!(p.z() < domain.zmin()) && p.z() < domain.zmax())
|
||||
return std::make_pair(p, Offset());
|
||||
|
||||
typename Geom_traits::Construct_point_3 cp = gt.construct_point_3_object();
|
||||
|
||||
// Numerical approximations might create inconsistencies between the constructions
|
||||
// and the comparisons. For example in a cubic domain of size 2:
|
||||
// 1. initial point: P(2+1e-17, 0, 0)
|
||||
// 2. the function computes an offset(1, 0, 0),
|
||||
// 3. P + (-1, 0, 0) * domain_size constructs Q(-1e-17, 0, 0) // numerical approximation
|
||||
// 4. the function computes an offset of (-1, 0, 0)
|
||||
// 5. Q + (1, 0, 0) * domain_size constructs (2+1e-17, 0, 0) (that is P)
|
||||
// And the function is looping...
|
||||
//
|
||||
// If this is happening the 'Last_change' enum will break this infinite
|
||||
// loop and return the wrong point and the 'encountered_issue' bool will be
|
||||
// set to 'true'. An exact version of this function should then be called.
|
||||
|
||||
enum Last_change {
|
||||
NO_LAST_CHANGE,
|
||||
INCREASED_X, DECREASED_X, INCREASED_Y, DECREASED_Y, INCREASED_Z, DECREASED_Z
|
||||
};
|
||||
|
||||
Last_change lc = NO_LAST_CHANGE;
|
||||
bool in = false;
|
||||
|
||||
Offset transl(0, 0, 0);
|
||||
Point dp;
|
||||
|
||||
while(!in)
|
||||
{
|
||||
dp = cp(p, transl);
|
||||
|
||||
if(dp.x() < domain.xmin())
|
||||
{
|
||||
if(lc == DECREASED_X) // stuck in a loop
|
||||
break;
|
||||
|
||||
lc = INCREASED_X;
|
||||
transl.x() += 1;
|
||||
}
|
||||
else if(dp.y() < domain.ymin())
|
||||
{
|
||||
if(lc == DECREASED_Y) // stuck in a loop
|
||||
break;
|
||||
|
||||
lc = INCREASED_Y;
|
||||
transl.y() += 1;
|
||||
}
|
||||
else if(dp.z() < domain.zmin())
|
||||
{
|
||||
if(lc == DECREASED_Z) // stuck in a loop
|
||||
break;
|
||||
|
||||
lc = INCREASED_Z;
|
||||
transl.z() += 1;
|
||||
}
|
||||
else if(!(dp.x() < domain.xmax()))
|
||||
{
|
||||
if(lc == INCREASED_X) // stuck in a loop
|
||||
break;
|
||||
|
||||
lc = DECREASED_X;
|
||||
transl.x() -= 1;
|
||||
}
|
||||
else if(!(dp.y() < domain.ymax()))
|
||||
{
|
||||
if(lc == INCREASED_Y) // stuck in a loop
|
||||
break;
|
||||
|
||||
lc = DECREASED_Y;
|
||||
transl.y() -= 1;
|
||||
}
|
||||
else if(!(dp.z() < domain.zmax()))
|
||||
{
|
||||
if(lc == INCREASED_Z) // stuck in a loop
|
||||
break;
|
||||
|
||||
lc = DECREASED_Z;
|
||||
transl.z() -= 1;
|
||||
}
|
||||
else
|
||||
{
|
||||
in = true;
|
||||
}
|
||||
}
|
||||
|
||||
std::pair<Point, Offset> pp(p, transl);
|
||||
|
||||
if(dp.x() < domain.xmin() || !(dp.x() < domain.xmax()) ||
|
||||
dp.y() < domain.ymin() || !(dp.y() < domain.ymax()) ||
|
||||
dp.z() < domain.zmin() || !(dp.z() < domain.zmax()))
|
||||
{
|
||||
encountered_issue = true;
|
||||
pp = construct_periodic_point_exact(p, gt);
|
||||
}
|
||||
|
||||
return pp;
|
||||
}
|
||||
|
||||
template <typename Gt_>
|
||||
bool
|
||||
is_point_too_close_to_border(const std::pair<typename Gt_::Point_3,
|
||||
typename Gt_::Periodic_3_offset_3>& pbp,
|
||||
const Gt_& gt)
|
||||
{
|
||||
typedef Gt_ Geom_traits;
|
||||
typedef typename Geom_traits::FT FT;
|
||||
typedef typename Geom_traits::Point_3 Bare_point;
|
||||
typedef typename Geom_traits::Iso_cuboid_3 Iso_cuboid;
|
||||
|
||||
typename Geom_traits::Construct_point_3 cp = gt.construct_point_3_object();
|
||||
|
||||
const Bare_point p = cp(pbp.first /*point*/, pbp.second /*offset*/);
|
||||
const FT px = p.x();
|
||||
const FT py = p.y();
|
||||
const FT pz = p.z();
|
||||
|
||||
const Iso_cuboid& domain = gt.get_domain();
|
||||
const FT dxm = domain.xmin();
|
||||
const FT dym = domain.ymin();
|
||||
const FT dzm = domain.zmin();
|
||||
const FT dxM = domain.xmax();
|
||||
const FT dyM = domain.ymax();
|
||||
const FT dzM = domain.zmax();
|
||||
|
||||
// simply comparing to FT::epsilon() is probably not completely satisfactory
|
||||
const FT eps = std::numeric_limits<FT>::epsilon();
|
||||
|
||||
FT diff = CGAL::abs(px - dxm);
|
||||
if(diff < eps && diff > 0) return true;
|
||||
diff = CGAL::abs(px - dxM);
|
||||
if(diff < eps && diff > 0) return true;
|
||||
diff = CGAL::abs(py - dym);
|
||||
if(diff < eps && diff > 0) return true;
|
||||
diff = CGAL::abs(py - dyM);
|
||||
if(diff < eps && diff > 0) return true;
|
||||
diff = CGAL::abs(pz - dzm);
|
||||
if(diff < eps && diff > 0) return true;
|
||||
diff = CGAL::abs(pz - dzM);
|
||||
if(diff < eps && diff > 0) return true;
|
||||
return false;
|
||||
}
|
||||
|
||||
template <typename Gt_>
|
||||
typename Gt_::Point_3
|
||||
snap_to_domain_border(const typename Gt_::Point_3& p, const Gt_& gt)
|
||||
{
|
||||
typedef Gt_ Geom_traits;
|
||||
typedef typename Geom_traits::FT FT;
|
||||
typedef typename Geom_traits::Iso_cuboid_3 Iso_cuboid;
|
||||
|
||||
const FT px = p.x();
|
||||
const FT py = p.y();
|
||||
const FT pz = p.z();
|
||||
FT sx = px, sy = py, sz = pz;
|
||||
|
||||
const Iso_cuboid& domain = gt.get_domain();
|
||||
const FT dxm = domain.xmin();
|
||||
const FT dym = domain.ymin();
|
||||
const FT dzm = domain.zmin();
|
||||
const FT dxM = domain.xmax();
|
||||
const FT dyM = domain.ymax();
|
||||
const FT dzM = domain.zmax();
|
||||
|
||||
// simply comparing to FT::epsilon() is probably not completely satisfactory
|
||||
const FT eps = std::numeric_limits<FT>::epsilon();
|
||||
|
||||
if(CGAL::abs(px - dxm) < eps) sx = dxm;
|
||||
if(CGAL::abs(px - dxM) < eps) sx = dxM;
|
||||
if(CGAL::abs(py - dym) < eps) sy = dym;
|
||||
if(CGAL::abs(py - dyM) < eps) sy = dyM;
|
||||
if(CGAL::abs(pz - dzm) < eps) sz = dzm;
|
||||
if(CGAL::abs(pz - dzM) < eps) sz = dzM;
|
||||
|
||||
std::cout << "epsilon: " << eps << std::endl;
|
||||
std::cout.precision(20);
|
||||
std::cout << "snapped " << p << " to " << sx << " " << sy << " " << sz << std::endl;
|
||||
|
||||
return gt.construct_point_3_object()(sx, sy, sz);
|
||||
}
|
||||
|
||||
template <typename Gt_>
|
||||
typename Gt_::Weighted_point_3
|
||||
snap_to_domain_border(const typename Gt_::Weighted_point_3& p,
|
||||
const Gt_& gt)
|
||||
{
|
||||
typedef Gt_ Geom_traits;
|
||||
typedef typename Geom_traits::Point_3 Bare_point;
|
||||
|
||||
typename Geom_traits::Compute_weight_3 cw = gt.compute_weight_3_object();
|
||||
|
||||
const Bare_point snapped_p = snap_to_domain_border(gt.construct_point_3_object()(p), gt);
|
||||
|
||||
return gt.construct_weighted_point_3_object()(snapped_p, cw(p));
|
||||
}
|
||||
|
||||
/// transform a bare point (living anywhere in space) into the canonical
|
||||
/// instance of the same bare point that lives inside the base domain
|
||||
template <typename Gt_>
|
||||
typename Gt_::Point_3
|
||||
robust_canonicalize_point(const typename Gt_::Point_3& p, const Gt_& gt)
|
||||
{
|
||||
typedef Gt_ Geom_traits;
|
||||
typedef typename Geom_traits::Point_3 Bare_point;
|
||||
typedef typename Geom_traits::Weighted_point_3 Weighted_point;
|
||||
typedef typename Geom_traits::Periodic_3_offset_3 Offset;
|
||||
typedef typename Geom_traits::Iso_cuboid_3 Iso_cuboid;
|
||||
|
||||
const Iso_cuboid& domain = gt.get_domain();
|
||||
if(p.x() >= domain.xmin() && p.x() < domain.xmax() &&
|
||||
p.y() >= domain.ymin() && p.y() < domain.ymax() &&
|
||||
p.z() >= domain.zmin() && p.z() < domain.zmax())
|
||||
return p;
|
||||
|
||||
bool should_snap = false;
|
||||
std::pair<Bare_point, Offset> pbp = construct_periodic_point(p, should_snap, gt);
|
||||
|
||||
if(!should_snap)
|
||||
{
|
||||
// Even if there is no issue while constructing the canonical point,
|
||||
// snap the point if it's too close to a border of the domain
|
||||
should_snap = is_point_too_close_to_border(pbp, gt);
|
||||
}
|
||||
|
||||
if(should_snap)
|
||||
{
|
||||
Bare_point sp = snap_to_domain_border(p, gt);
|
||||
|
||||
// might have snapped to a 'max' of the domain, which is not in the domain
|
||||
// note: we could snap to 'min' all the time in 'snap_to_domain_border'
|
||||
// but this is clearer like that (and costs very little since we should
|
||||
// not have to use exact computations too often)
|
||||
return robust_canonicalize_point(sp, gt);
|
||||
}
|
||||
|
||||
typename Geom_traits::Construct_point_3 cp = gt.construct_point_3_object();
|
||||
|
||||
Bare_point canonical_p = cp(pbp.first /*point*/, pbp.second /*offset*/);
|
||||
CGAL_postcondition( !(canonical_p.x() < domain.xmin()) &&
|
||||
(canonical_p.x() < domain.xmax()) );
|
||||
CGAL_postcondition( !(canonical_p.y() < domain.ymin()) &&
|
||||
(canonical_p.y() < domain.ymax()) );
|
||||
CGAL_postcondition( !(canonical_p.z() < domain.zmin()) &&
|
||||
(canonical_p.z() < domain.zmax()) );
|
||||
|
||||
return canonical_p;
|
||||
}
|
||||
|
||||
/// transform a weighted point (living anywhere in space) into the canonical
|
||||
/// instance of the same weighted point that lives inside the base domain
|
||||
template <typename Gt_>
|
||||
typename Gt_::Weighted_point_3
|
||||
robust_canonicalize_point(const typename Gt_::Weighted_point_3& wp, const Gt_& gt)
|
||||
{
|
||||
typedef Gt_ Geom_traits;
|
||||
typedef typename Geom_traits::Point_3 Bare_point;
|
||||
typedef typename Geom_traits::Weighted_point_3 Weighted_point;
|
||||
typedef typename Geom_traits::Iso_cuboid_3 Iso_cuboid;
|
||||
|
||||
const Iso_cuboid& domain = gt.get_domain();
|
||||
if(wp.x() >= domain.xmin() && wp.x() < domain.xmax() &&
|
||||
wp.y() >= domain.ymin() && wp.y() < domain.ymax() &&
|
||||
wp.z() >= domain.zmin() && wp.z() < domain.zmax())
|
||||
return wp;
|
||||
|
||||
typename Geom_traits::Construct_point_3 cp = gt.construct_point_3_object();
|
||||
typename Geom_traits::Compute_weight_3 cw = gt.compute_weight_3_object();
|
||||
typename Geom_traits::Construct_weighted_point_3 cwp = gt.construct_weighted_point_3_object();
|
||||
|
||||
const Bare_point& bp = cp(wp);
|
||||
Bare_point canonical_point = robust_canonicalize_point(bp, gt);
|
||||
|
||||
return cwp(canonical_point, cw(wp));
|
||||
}
|
||||
|
||||
} // end namespace internal
|
||||
} // end namespace Periodic_3_triangulation_3
|
||||
} // end namespace CGAL
|
||||
|
||||
#endif // CGAL_PERIODIC_3_TRIANGULATION_3_CANONICALIZE_HELPER_H
|
||||
Loading…
Reference in New Issue