mirror of https://github.com/CGAL/cgal
adding predicates and filtering
This commit is contained in:
Olivier Devillers
parent
76e95c9689
commit
0086efbf49
|
|
@ -1716,7 +1716,12 @@ Kernel_23/include/CGAL/internal/Projection_traits_3.h -text
|
|||
Kernel_23/test/Kernel_23/CMakeLists.txt -text
|
||||
Kernel_d/doc_tex/Kernel_d/hypercube.png -text
|
||||
Kernel_d/doc_tex/Kernel_d_ref/Kernel_Compute_coordinate_d.tex -text
|
||||
Kernel_d/include/CGAL/Filtered_kernel_d.h -text
|
||||
Kernel_d/include/CGAL/Kernel_d/Cartesian_const_iterator_d.h -text
|
||||
Kernel_d/include/CGAL/Kernel_d/Cartesian_converter_d.h -text
|
||||
Kernel_d/include/CGAL/Kernel_d/Interval_linear_algebra.h -text
|
||||
Kernel_d/include/CGAL/Kernel_d/interface_macros_d.h -text
|
||||
Kernel_d/include/CGAL/Referenced_argument.h -text
|
||||
Kernel_d/test/Kernel_d/Linear_algebra-test.cmd eol=lf
|
||||
Kinetic_data_structures/demo/Kinetic_data_structures/data/after002 -text
|
||||
Kinetic_data_structures/demo/Kinetic_data_structures/data/after010 -text
|
||||
|
|
|
|||
|
|
@ -242,7 +242,9 @@ public:
|
|||
typedef Position_on_lineCd<Self> Position_on_line_d;
|
||||
typedef Barycentric_coordinatesCd<Self> Barycentric_coordinates_d;
|
||||
typedef OrientationCd<Self> Orientation_d;
|
||||
typedef Coaffine_orientationCd<Self> Coaffine_orientation_d;
|
||||
typedef Side_of_oriented_sphereCd<Self> Side_of_oriented_sphere_d;
|
||||
typedef Side_of_oriented_subsphereCd<Self> Side_of_oriented_subsphere_d;
|
||||
typedef Side_of_bounded_sphereCd<Self> Side_of_bounded_sphere_d;
|
||||
typedef Contained_in_simplexCd<Self> Contained_in_simplex_d;
|
||||
typedef Contained_in_affine_hullCd<Self> Contained_in_affine_hull_d;
|
||||
|
|
@ -274,8 +276,12 @@ public:
|
|||
{ return Barycentric_coordinates_d(); }
|
||||
Orientation_d orientation_d_object() const
|
||||
{ return Orientation_d(); }
|
||||
Coaffine_orientation_d coaffine_orientation_d_object() const
|
||||
{ return Coaffine_orientation_d(); }
|
||||
Side_of_oriented_sphere_d side_of_oriented_sphere_d_object() const
|
||||
{ return Side_of_oriented_sphere_d(); }
|
||||
Side_of_oriented_subsphere_d side_of_oriented_subsphere_d_object() const
|
||||
{ return Side_of_oriented_subsphere_d(); }
|
||||
Side_of_bounded_sphere_d side_of_bounded_sphere_d_object() const
|
||||
{ return Side_of_bounded_sphere_d(); }
|
||||
Contained_in_simplex_d contained_in_simplex_d_object() const
|
||||
|
|
|
|||
|
|
@ -0,0 +1,54 @@
|
|||
|
||||
#ifndef CGAL_FILTERED_KERNEL_D_H
|
||||
#define CGAL_FILTERED_KERNEL_D_H
|
||||
|
||||
#include <CGAL/Filtered_predicate.h>
|
||||
#include <CGAL/internal/Exact_type_selector.h>
|
||||
#include <CGAL/Kernel_d/Cartesian_converter_d.h>
|
||||
|
||||
#include<boost/pool/pool_alloc.hpp>
|
||||
|
||||
namespace CGAL {
|
||||
|
||||
template<typename Kernel> // a dD kernel we want to filter
|
||||
struct Filtered_kernel_d : public Cartesian_d<typename Kernel::FT>
|
||||
{
|
||||
typedef typename Kernel::LA LA;
|
||||
typedef typename Kernel::RT RT; // Ring type
|
||||
typedef typename Kernel::FT FT; // Field type
|
||||
|
||||
typedef Cartesian_d<FT> Base;
|
||||
typedef Filtered_kernel_d<Kernel> Self;
|
||||
|
||||
// an exact number type
|
||||
typedef typename internal::Exact_type_selector<RT>::Type Exact_nt;
|
||||
|
||||
// the corresponding exact kernel
|
||||
//typedef Linear_algebraCd< Exact_nt, boost::pool_allocator<Exact_nt> > Exact_linalg;
|
||||
typedef Linear_algebraCd< Exact_nt > Exact_linalg;
|
||||
typedef Cartesian_d<Exact_nt, Exact_linalg> EK;
|
||||
|
||||
// the kernel used for filtered predicates
|
||||
typedef Interval_nt<false> IA;
|
||||
//typedef Linear_algebraCd<IA, boost::pool_allocator<IA> > Interval_linalg;
|
||||
typedef Linear_algebraCd<IA > Interval_linalg;
|
||||
typedef Cartesian_d<IA, Interval_linalg > FK;
|
||||
|
||||
// the converter
|
||||
typedef Cartesian_converter_d<Base, EK> C2E;
|
||||
typedef Cartesian_converter_d<Base, FK> C2F;
|
||||
|
||||
// we change the predicates.
|
||||
#define CGAL_Kernel_pred(P, Pf) \
|
||||
typedef Filtered_predicate<typename EK::P, typename FK::P, C2E, C2F> P; \
|
||||
P Pf() const { return P(); }
|
||||
|
||||
// we don't touch the constructions.
|
||||
#define CGAL_Kernel_cons(Y,Z)
|
||||
|
||||
#include <CGAL/Kernel_d/interface_macros_d.h>
|
||||
};
|
||||
|
||||
} //namespace CGAL
|
||||
|
||||
#endif // CGAL_FILTERED_KERNEL_D_H
|
||||
|
|
@ -0,0 +1,261 @@
|
|||
// Copyright (c) 2001-2004 Utrecht University (The Netherlands),
|
||||
// ETH Zurich (Switzerland), Freie Universitaet Berlin (Germany),
|
||||
// INRIA Sophia-Antipolis (France), Martin-Luther-University Halle-Wittenberg
|
||||
// (Germany), Max-Planck-Institute Saarbruecken (Germany), RISC Linz (Austria),
|
||||
// and Tel-Aviv University (Israel). 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 Lesser General Public License as
|
||||
// published by the Free Software Foundation; version 2.1 of the License.
|
||||
// See the file LICENSE.LGPL distributed with CGAL.
|
||||
//
|
||||
// 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) : Sylvain Pion <Sylvain.Pion@sophia.inria.fr>
|
||||
// Menelaos Karavelas <mkaravel@cse.nd.edu>
|
||||
// Samuel Hornus <Samuel.Hornus@sophia.inria.fr>
|
||||
|
||||
#ifndef CGAL_CARTESIAN_CONVERTER_D_H
|
||||
#define CGAL_CARTESIAN_CONVERTER_D_H
|
||||
|
||||
// This file contains the definition of a dD kernel converter, based on Cartesian
|
||||
// representation. It should work between *Cartesian_d<A> and *Cartesian_d<B>,
|
||||
// provided you give a NT converter from A to B.
|
||||
// TODO: There's NO Homogeneous counterpart at the moment.
|
||||
// TODO: Maybe we can derive Cartesian_converter_d from Cartesian_converter ?
|
||||
// But I don't know if it is interesting or not.
|
||||
|
||||
#include <CGAL/basic.h>
|
||||
#include <CGAL/NT_converter.h>
|
||||
#include <CGAL/Enum_converter.h>
|
||||
#include <CGAL/Iterator_transform.h>
|
||||
#include <vector>
|
||||
|
||||
#include <CGAL/Referenced_argument.h>
|
||||
|
||||
namespace CGAL {
|
||||
|
||||
// Guess which compiler needs this work around ?
|
||||
namespace CGALi {
|
||||
template < typename K1, typename K2 >
|
||||
struct Default_converter_d {
|
||||
typedef typename K1::FT FT1;
|
||||
typedef typename K2::FT FT2;
|
||||
typedef ::CGAL::NT_converter<FT1, FT2> Type;
|
||||
};
|
||||
} // namespace CGALi
|
||||
|
||||
template < class K1, class K2,
|
||||
//class Converter = NT_converter<typename K1::FT, typename K2::FT> >
|
||||
class Converter = typename CGALi::Default_converter_d<K1, K2>::Type >
|
||||
class Cartesian_converter_d : public Enum_converter
|
||||
{
|
||||
typedef Enum_converter Base;
|
||||
typedef Cartesian_converter_d<K1, K2, Converter> Self;
|
||||
|
||||
public:
|
||||
typedef K1 Source_kernel;
|
||||
typedef K2 Target_kernel;
|
||||
typedef Converter Number_type_converter;
|
||||
|
||||
#ifdef CGAL_CFG_USING_BASE_MEMBER_BUG
|
||||
bool operator()(bool b) const { return Base::operator()(b); }
|
||||
Sign operator()(Sign s) const { return Base::operator()(s); }
|
||||
|
||||
Oriented_side operator()(Oriented_side os) const {
|
||||
return Base::operator()(os);
|
||||
}
|
||||
|
||||
Bounded_side operator()(Bounded_side bs) const {
|
||||
return Base::operator()(bs);
|
||||
}
|
||||
|
||||
Comparison_result operator()(Comparison_result cr) const {
|
||||
return Base::operator()(cr);
|
||||
}
|
||||
|
||||
Angle operator()(Angle a) const { return Base::operator()(a); }
|
||||
#else
|
||||
using Base::operator();
|
||||
#endif
|
||||
|
||||
Cartesian_converter_d() // To shut up a warning with SunPRO.
|
||||
: c(), k(), result_point_(20) {}
|
||||
|
||||
Origin
|
||||
operator()(const Origin& o) const
|
||||
{
|
||||
return o;
|
||||
}
|
||||
|
||||
Null_vector
|
||||
operator()(const Null_vector& n) const
|
||||
{
|
||||
return n;
|
||||
}
|
||||
|
||||
typename K2::FT
|
||||
operator()(const typename K1::FT &a) const
|
||||
{
|
||||
return c(a);
|
||||
}
|
||||
|
||||
std::vector<Object>
|
||||
operator()(const std::vector<Object>& v) const
|
||||
{
|
||||
std::vector<Object> res;
|
||||
res.reserve(v.size());
|
||||
for(unsigned int i = 0; i < v.size(); i++) {
|
||||
res.push_back(operator()(v[i]));
|
||||
}
|
||||
return res;
|
||||
}
|
||||
|
||||
int operator()(const int &a)
|
||||
{
|
||||
return a;
|
||||
}
|
||||
|
||||
typename K2::Point_d
|
||||
operator()(const typename K1::Point_d &a) const
|
||||
{
|
||||
typedef typename K2::Point_d Point_d;
|
||||
typedef typename K1::Point_d::Cartesian_const_iterator Coord_iter;
|
||||
typedef Iterator_transform<Coord_iter, Converter> It;
|
||||
return Point_d(a.dimension(), It(a.cartesian_begin()), It(a.cartesian_end()));
|
||||
}
|
||||
|
||||
typename K2::Vector_d
|
||||
operator()(const typename K1::Vector_d &a) const
|
||||
{
|
||||
typedef typename K2::Vector_d Vector_d;
|
||||
typedef typename K1::Point_d::Cartesian_const_iterator Coord_iter;
|
||||
typedef Iterator_transform<Coord_iter, Converter> It;
|
||||
return Vector_d(a.dimension(), It(a.cartesian_begin()), It(a.cartesian_end()));
|
||||
}
|
||||
|
||||
typename K2::Direction_d
|
||||
operator()(const typename K1::Direction_d &a) const
|
||||
{
|
||||
typedef typename K2::Direction_d Direction_d;
|
||||
return Direction_d(operator()(a.vector()));
|
||||
}
|
||||
|
||||
typename K2::Segment_d
|
||||
operator()(const typename K1::Segment_d &a) const
|
||||
{
|
||||
typedef typename K2::Segment_d Segment_d;
|
||||
return Segment_d(operator()(a.source()), operator()(a.target()));
|
||||
}
|
||||
|
||||
typename K2::Line_d
|
||||
operator()(const typename K1::Line_d &a) const
|
||||
{
|
||||
typedef typename K2::Line_d Line_d;
|
||||
return Line_d(operator()(a.point(0)), operator()(a.direction()));
|
||||
}
|
||||
|
||||
typename K2::Ray_d
|
||||
operator()(const typename K1::Ray_d &a) const
|
||||
{
|
||||
typedef typename K2::Ray_d Ray_d;
|
||||
return Ray_d(operator()(a.source()), operator()(a.second_point()));
|
||||
}
|
||||
|
||||
typename K2::Sphere_d
|
||||
operator()(const typename K1::Sphere_d &a) const
|
||||
{
|
||||
typedef typename K2::Sphere_d Sphere_d;
|
||||
typedef typename K1::Sphere_d::point_iterator Coord_iter;
|
||||
// TODO: Check that the use of Iterator_transform is correct
|
||||
typedef Iterator_transform<Coord_iter, Converter> It;
|
||||
return Sphere_d(a.dimension(), It(a.points_begin()), It(a.points_end()));
|
||||
}
|
||||
|
||||
typename K2::Hyperplane_d
|
||||
operator()(const typename K1::Hyperplane_d &a) const
|
||||
{
|
||||
typedef typename K2::Hyperplane_d Hyperplane_d;
|
||||
typedef typename K1::Hyperplane_d::Coefficient_const_iterator Coord_iter;
|
||||
// TODO: Check that the use of Iterator_transform is correct
|
||||
typedef Iterator_transform<Coord_iter, Converter> It;
|
||||
return Hyperplane_d(a.dimension(), It(a.coefficients_begin()), It(a.coefficients_end()));
|
||||
}
|
||||
|
||||
typename K2::Iso_box_d
|
||||
operator()(const typename K1::Iso_box_d &a) const
|
||||
{
|
||||
typedef typename K2::Iso_box_d Iso_box_d;
|
||||
return Iso_box_d(operator()((a.min)()), operator()((a.max)()));
|
||||
}
|
||||
|
||||
/*std::vector<int> &
|
||||
operator()(const std::vector<int> & a) const
|
||||
{
|
||||
return const_cast<std::vector<int> &>(a);
|
||||
}*/
|
||||
std::vector<int>::iterator
|
||||
operator()(const std::vector<int>::iterator & a) const
|
||||
{
|
||||
return a;
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
Referenced_argument<T> & operator()(const Referenced_argument<T> & a) const
|
||||
{
|
||||
return const_cast<Referenced_argument<T> &> (a);
|
||||
}
|
||||
|
||||
// ---- The following allows to convert iterators to Point_d, Vector_d,
|
||||
// etc... by adapting the iterator with an 'on-the-fly' converter. This is
|
||||
// not the most efficient because, if the iterator is dereferenced N times,
|
||||
// its value will be converted N times... TODO: Convert once and for all:
|
||||
// this should probably be done in Filtered_kernel.h (or a new
|
||||
// Filtered_kernel_d.h file).
|
||||
template<typename IterBase>
|
||||
struct Specialized_converter : public Self
|
||||
{
|
||||
typedef typename IterBase::value_type argument_type;
|
||||
typedef typename argument_type::
|
||||
template WithAnotherKernel<K2>::Type result_type;
|
||||
};
|
||||
|
||||
template<typename IterBase>
|
||||
Iterator_transform<IterBase, Specialized_converter<IterBase> >
|
||||
operator()(const IterBase & a) const
|
||||
{
|
||||
return Iterator_transform<IterBase, Specialized_converter<IterBase> >(a);
|
||||
}
|
||||
// ---- ---- ---- ---- ---- ---- //
|
||||
|
||||
private:
|
||||
Converter c;
|
||||
K2 k;
|
||||
typename K2::Point_d result_point_;
|
||||
};
|
||||
|
||||
// Specialization when converting to the same kernel,
|
||||
// to avoid making copies.
|
||||
template < class K, class C >
|
||||
class Cartesian_converter_d <K, K, C>
|
||||
{
|
||||
public:
|
||||
typedef K Source_kernel;
|
||||
typedef K Target_kernel;
|
||||
typedef C Number_type_converter;
|
||||
|
||||
template < typename T >
|
||||
const T& operator()(const T&t) const { return t; }
|
||||
};
|
||||
|
||||
} //namespace CGAL
|
||||
|
||||
#endif // CGAL_CARTESIAN_CONVERTER_D_H
|
||||
|
|
@ -0,0 +1,94 @@
|
|||
|
||||
#ifndef CGAL_INTERVAL_LINEAR_ALGEBRA_H
|
||||
#define CGAL_INTERVAL_LINEAR_ALGEBRA_H
|
||||
|
||||
#include <CGAL/Interval_nt.h>
|
||||
#include <CGAL/predicates/sign_of_determinant.h>
|
||||
|
||||
namespace CGAL {
|
||||
|
||||
namespace CGALi
|
||||
{
|
||||
|
||||
template<class Matrix>
|
||||
Sign // The matrix is row major: M[i] represents row i.
|
||||
sign_of_determinantDxD_with_interval_arithmetic(Matrix & M)
|
||||
// attempts to compute the determinant using interval arithmetic
|
||||
{
|
||||
typedef Interval_nt<false> NT; // Field number type
|
||||
int sign = 1;
|
||||
int curdim = M.dimension().first;
|
||||
for(int col = 0; col < curdim; ++col)
|
||||
{
|
||||
int pivot_line = -1;
|
||||
NT pivot = 0.0;
|
||||
for(int l = col; l < curdim; ++l)
|
||||
{
|
||||
NT candidate = CGAL::abs(M[l][col]);
|
||||
if( CGAL::certainly(0.0 < candidate) )
|
||||
if( pivot.inf() < candidate.inf() )
|
||||
{
|
||||
pivot_line = l;
|
||||
pivot = candidate;
|
||||
}
|
||||
}
|
||||
if( -1 == pivot_line )
|
||||
{
|
||||
throw CGAL::Uncertain_conversion_exception("undecidable interval computation of determinant");
|
||||
}
|
||||
// if the pivot value is negative, invert the sign
|
||||
if( M[pivot_line][col] < 0.0 )
|
||||
sign = - sign;
|
||||
// put the pivot cell on the diagonal
|
||||
if( pivot_line > col )
|
||||
{
|
||||
std::swap(M[col], M[pivot_line]); // this takes constant time with std::vector<>
|
||||
sign = - sign;
|
||||
}
|
||||
// using matrix-line operations, set zero in all cells below the pivot cell.
|
||||
// This costs roughly (curdim-col-1)^2 mults and adds, because we don't actually
|
||||
// compute the zeros below the pivot cell
|
||||
NT denom = NT(1.0) / M[col][col];
|
||||
for(int line = col + 1; line < curdim; ++line)
|
||||
{
|
||||
NT numer = M[line][col] * denom;
|
||||
for (int c = col + 1; c < curdim; ++c)
|
||||
M[line][c] -= numer * M[col][c];
|
||||
}
|
||||
}
|
||||
return Sign(sign);
|
||||
}
|
||||
|
||||
} // end of namespace CGALi
|
||||
|
||||
template<>
|
||||
inline
|
||||
Sign
|
||||
Linear_algebraCd<Interval_nt_advanced>::sign_of_determinant(const Matrix & M)
|
||||
{
|
||||
switch( M.dimension().first )
|
||||
{
|
||||
case 2:
|
||||
return CGAL::sign_of_determinant( M(0,0), M(0,1),
|
||||
M(1,0), M(1,1));
|
||||
break;
|
||||
case 3:
|
||||
return CGAL::sign_of_determinant( M(0,0), M(0,1), M(0,2),
|
||||
M(1,0), M(1,1), M(1,2),
|
||||
M(2,0), M(2,1), M(2,2));
|
||||
break;
|
||||
case 4:
|
||||
return CGAL::sign_of_determinant( M(0,0), M(0,1), M(0,2), M(0,3),
|
||||
M(1,0), M(1,1), M(1,2), M(1,3),
|
||||
M(2,0), M(2,1), M(2,2), M(2,3),
|
||||
M(3,0), M(3,1), M(3,2), M(3,3));
|
||||
break;
|
||||
default:
|
||||
return CGALi::sign_of_determinantDxD_with_interval_arithmetic(const_cast<Matrix &>(M));
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
} //namespace CGAL
|
||||
|
||||
#endif // CGAL_INTERVAL_LINEAR_ALGEBRA_H
|
||||
|
|
@ -26,6 +26,7 @@
|
|||
|
||||
#include <CGAL/basic.h>
|
||||
#include <CGAL/enum.h>
|
||||
#include <CGAL/Referenced_argument.h>
|
||||
|
||||
#undef CGAL_KD_TRACE
|
||||
#undef CGAL_KD_TRACEN
|
||||
|
|
@ -212,6 +213,91 @@ Orientation operator()(ForwardIterator first, ForwardIterator last)
|
|||
}
|
||||
};
|
||||
|
||||
/* This predicates tests the orientation of (k+1) points that span a
|
||||
* k-dimensional affine subspace of the ambiant d-dimensional space. We
|
||||
* greedily search for an orthogonal projection on a k-dim axis aligned
|
||||
* subspace on which the (full k-dim) predicates answers POSITIVE or NEGATIVE.
|
||||
* If no such subspace is found, return COPLANAR.
|
||||
* IMPORTANT TODO: Current implementation is VERY bad with filters: if one
|
||||
* determinant fails in the filtering step, then all the subsequent ones wil be
|
||||
* in exact arithmetic :-(
|
||||
* TODO: store the axis-aligned subspace that was found in order to avoid
|
||||
* re-searching for it for subsequent calls to operator()
|
||||
*/
|
||||
template <class R>
|
||||
struct Coaffine_orientationCd
|
||||
{
|
||||
typedef typename R::Point_d Point_d;
|
||||
typedef typename R::LA LA;
|
||||
typedef Orientation result_type;
|
||||
// typedef internal::stateful_predicate_tag predicate_category;
|
||||
typedef std::vector<int> Axes;
|
||||
struct State
|
||||
{
|
||||
Axes axes_;
|
||||
bool axes_found_;
|
||||
State(bool b) : axes_(), axes_found_(b) {}
|
||||
};
|
||||
mutable State state_;
|
||||
//typedef Referenced_argument<std::vector<int> > Axes;
|
||||
//typedef Referenced_argument<bool> Ref_bool;
|
||||
|
||||
Coaffine_orientationCd() : state_(false) {}
|
||||
State & state() { return state_; }
|
||||
const State & state() const { return state_; }
|
||||
|
||||
template <class ForwardIterator>
|
||||
result_type operator()(ForwardIterator first, ForwardIterator last) const
|
||||
{
|
||||
TUPLE_DIM_CHECK(first,last,Coaffine_orientation_d);
|
||||
// |k| is the dimension of the affine subspace
|
||||
const int k = std::distance(first,last) - 1;
|
||||
// |d| is the dimension of the ambiant space
|
||||
int d = first->dimension();
|
||||
CGAL_assertion_msg(k <= d,
|
||||
"Coaffine_orientation_d: needs less that (first->dimension() + 1) points.");
|
||||
if( false == state_.axes_found_ )
|
||||
{
|
||||
state_.axes_.resize(d + 1);
|
||||
// We start by choosing the first |k| axes to define a plane of projection
|
||||
int i = 0;
|
||||
for(; i < k; ++i) state_.axes_[i] = i;
|
||||
for(; i < d + 1; ++i) state_.axes_[i] = -1;
|
||||
}
|
||||
typename ForwardIterator::value_type l = *(first + k);
|
||||
//ForwardIterator l = first + k;
|
||||
typename LA::Matrix M(k); // quadratic
|
||||
while( true )
|
||||
{
|
||||
for (int i = 0; i < k; ++i)
|
||||
{
|
||||
typename ForwardIterator::value_type fpi = *(first + i);
|
||||
for (int j = 0; j < k; ++j)
|
||||
M(i,j) = fpi.cartesian(state_.axes_[j]) -
|
||||
l.cartesian(state_.axes_[j]);
|
||||
}
|
||||
Orientation o = Orientation(LA::sign_of_determinant(M));
|
||||
if( ( o != COPLANAR ) || state_.axes_found_ )
|
||||
{
|
||||
state_.axes_found_ = true;
|
||||
return o;
|
||||
}
|
||||
|
||||
// for generating all possible unordered k-uple in the range
|
||||
// [0 .. d-1]... we go to the next unordered k-uple
|
||||
int index = k - 1;
|
||||
while( (index >= 0) && (state_.axes_[index] == d - k + index) )
|
||||
--index;
|
||||
if( index < 0 )
|
||||
break;
|
||||
++state_.axes_[index];
|
||||
for(int i = 1; i < k - index; ++i)
|
||||
state_.axes_[index + i] = state_.axes_[index] + i;
|
||||
}
|
||||
return COPLANAR;
|
||||
}
|
||||
};
|
||||
|
||||
template <class R>
|
||||
struct Side_of_oriented_sphereCd {
|
||||
typedef typename R::Point_d Point_d;
|
||||
|
|
@ -247,6 +333,128 @@ Oriented_side operator()(ForwardIterator first, ForwardIterator last,
|
|||
}
|
||||
};
|
||||
|
||||
|
||||
/* This predicates takes k+1 points defining a k-sphere in d-dim space, and a
|
||||
* point |x| (assumed to lie in the same affine subspace spanned by the
|
||||
* k-sphere). It tests wether the point |x| lies in the positive or negative
|
||||
* side of the k-sphere.
|
||||
* The parameter |axis| contains the indices of k axis of the canonical base of
|
||||
* R^d, on which the affine subspace projects homeomorphically. We can thus
|
||||
* "complete" the k+1 points with d-k other points along the "non-used" axes
|
||||
* and then call the usual Side_of_oriented_sphereCd predicate.
|
||||
*/
|
||||
template < class R >
|
||||
struct Side_of_oriented_subsphereCd
|
||||
{
|
||||
typedef typename R::Point_d Point;
|
||||
typedef typename R::LA LA;
|
||||
typedef typename R::FT FT;
|
||||
typedef Oriented_side result_type;
|
||||
// typedef internal::stateless_predicate_tag predicate_category;
|
||||
|
||||
typedef typename R::Side_of_oriented_sphere_d Side_of_oriented_sphere;
|
||||
typedef typename R::Coaffine_orientation_d Coaffine_orientation;
|
||||
typedef typename Coaffine_orientation::Axes Axes;
|
||||
|
||||
// DATA MEMBERS
|
||||
mutable Coaffine_orientation ori_;
|
||||
mutable typename LA::Matrix M; // a square matrix of size (D+1)x(D+1)
|
||||
mutable unsigned int adjust_sign_;
|
||||
|
||||
Side_of_oriented_subsphereCd()
|
||||
: ori_(R().coaffine_orientation_d_object()), M(), adjust_sign_(0) { }
|
||||
|
||||
template < class ForwardIterator >
|
||||
result_type operator()(ForwardIterator first, ForwardIterator last, const Point & q) const
|
||||
{
|
||||
const int d = first->dimension();
|
||||
const int k = std::distance(first, last) - 1; // dimension of affine subspace
|
||||
CGAL_assertion_msg( k <= d, "too much points in range.");
|
||||
if( k == d )
|
||||
{
|
||||
Side_of_oriented_sphere sos;
|
||||
return sos(first, last, q); // perhaps slap user on the back of the head here?
|
||||
}
|
||||
if( M.row_dimension() < d+1 )
|
||||
M = typename LA::Matrix(d+1);
|
||||
if( ! ori_.state().axes_found_ )
|
||||
{
|
||||
Orientation o = ori_(first, last);
|
||||
// FIXME: why these two lines below? I would remove them [sam]
|
||||
if( COPLANAR == o )
|
||||
return ON_ORIENTED_BOUNDARY;
|
||||
CGAL_assertion( o == POSITIVE );
|
||||
// Now we can setup the fixed part of the matrix:
|
||||
int a(0);
|
||||
int j(k);
|
||||
typename Axes::iterator axis = ori_.state().axes_.begin();
|
||||
while( j < d )
|
||||
{
|
||||
while( a == *axis )
|
||||
{
|
||||
++a; ++axis;
|
||||
}
|
||||
adjust_sign_ = ( adjust_sign_ + j + a ) % 2;
|
||||
int i(0);
|
||||
for( ; i < a; ++i )
|
||||
M(i, j) = FT(0);
|
||||
M(i++, j) = FT(1); // i.e.: M(a, j) = 1
|
||||
for( ; i < d; ++i )
|
||||
M(i, j) = FT(0);
|
||||
++j;
|
||||
++a;
|
||||
}
|
||||
}
|
||||
typename ForwardIterator::value_type p1 = *first;
|
||||
FT SumFirst(0); // squared length of first subsphere point, seen as vector.
|
||||
for( int i = 0; i < d; ++i )
|
||||
{
|
||||
FT ci = p1.cartesian(i);
|
||||
SumFirst += ci * ci;
|
||||
}
|
||||
int j(0); // iterates overs columns/subsphere points
|
||||
++first;
|
||||
while( first != last )
|
||||
{
|
||||
typename ForwardIterator::value_type v = *first;
|
||||
FT Sum = FT(0);
|
||||
for( int i = 0; i < d; ++i )
|
||||
{
|
||||
FT ci = v.cartesian(i);
|
||||
M(i, j) = ci - p1.cartesian(i);
|
||||
Sum += ci * ci;
|
||||
}
|
||||
M(d, j) = Sum - SumFirst;
|
||||
++first;
|
||||
++j;
|
||||
}
|
||||
int a(0);
|
||||
typename Axes::iterator axis = ori_.state().axes_.begin();
|
||||
while( j < d )
|
||||
{
|
||||
while( a == *axis )
|
||||
{
|
||||
++a; ++axis;
|
||||
}
|
||||
M(d, j) = FT(1) + FT(2) * p1.cartesian(a);
|
||||
++j;
|
||||
++a;
|
||||
}
|
||||
FT Sum = FT(0);
|
||||
for( int i = 0; i < d; ++i )
|
||||
{
|
||||
FT ci = q.cartesian(i);
|
||||
M(i, d) = ci - p1.cartesian(i);
|
||||
Sum += ci * ci;
|
||||
}
|
||||
M(d, d) = Sum - SumFirst;
|
||||
if( 0 == ( adjust_sign_ % 2 ) )
|
||||
return Oriented_side( - LA::sign_of_determinant( M ) );
|
||||
else
|
||||
return Oriented_side( LA::sign_of_determinant( M ) );
|
||||
}
|
||||
};
|
||||
|
||||
template <class R>
|
||||
struct Side_of_bounded_sphereCd {
|
||||
typedef typename R::Point_d Point_d;
|
||||
|
|
|
|||
|
|
@ -0,0 +1,109 @@
|
|||
// Copyright (c) 2000-2004 Utrecht University (The Netherlands),
|
||||
// ETH Zurich (Switzerland), Freie Universitaet Berlin (Germany),
|
||||
// INRIA Sophia-Antipolis (France), Martin-Luther-University Halle-Wittenberg
|
||||
// (Germany), Max-Planck-Institute Saarbruecken (Germany), RISC Linz (Austria),
|
||||
// and Tel-Aviv University (Israel). 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 Lesser General Public License as
|
||||
// published by the Free Software Foundation; version 2.1 of the License.
|
||||
// See the file LICENSE.LGPL distributed with CGAL.
|
||||
//
|
||||
// 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: svn+ssh://scm.gforge.inria.fr/svn/cgal/branches/CGAL-3.3-branch/Kernel_23/include/CGAL/Kernel/interface_macros.h $
|
||||
// $Id: interface_macros.h 34893 2006-10-24 05:24:31Z spion $
|
||||
//
|
||||
//
|
||||
// Author(s) : Herve Bronnimann, Sylvain Pion, Susan Hert
|
||||
|
||||
// This file is intentionally not protected against re-inclusion.
|
||||
// It's aimed at being included from within a kernel traits class, this
|
||||
// way we share more code.
|
||||
|
||||
// It is the responsability of the including file to correctly set the 2
|
||||
// macros CGAL_Kernel_pred, CGAL_Kernel_cons and CGAL_Kernel_obj.
|
||||
// And they are #undefed at the end of this file.
|
||||
|
||||
#ifndef CGAL_Kernel_pred
|
||||
# define CGAL_Kernel_pred(X, Y)
|
||||
#endif
|
||||
|
||||
#ifndef CGAL_Kernel_cons
|
||||
# define CGAL_Kernel_cons(X, Y)
|
||||
#endif
|
||||
|
||||
#ifndef CGAL_Kernel_obj
|
||||
# define CGAL_Kernel_obj(X)
|
||||
#endif
|
||||
|
||||
CGAL_Kernel_obj(Point_d)
|
||||
CGAL_Kernel_obj(Vector_d)
|
||||
CGAL_Kernel_obj(Direction_d)
|
||||
CGAL_Kernel_obj(Hyperplane_d)
|
||||
CGAL_Kernel_obj(Sphere_d)
|
||||
CGAL_Kernel_obj(Iso_box_d)
|
||||
CGAL_Kernel_obj(Segment_d)
|
||||
CGAL_Kernel_obj(Ray_d)
|
||||
CGAL_Kernel_obj(Line_d)
|
||||
|
||||
|
||||
CGAL_Kernel_pred(Affinely_independent_d,
|
||||
affinely_independent_d_object)
|
||||
CGAL_Kernel_pred(Affine_rank_d,
|
||||
affine_rank_d_object)
|
||||
CGAL_Kernel_pred(Compare_lexicographically_d,
|
||||
compare_lexicographically_d_object)
|
||||
CGAL_Kernel_pred(Contained_in_affine_hull_d,
|
||||
contained_in_affine_hull_d_object)
|
||||
CGAL_Kernel_pred(Contained_in_linear_hull_d,
|
||||
contained_in_linear_hull_d_object)
|
||||
CGAL_Kernel_pred(Contained_in_simplex_d,
|
||||
contained_in_simplex_d_object)
|
||||
// TODO: create a Do_intersect_d functor
|
||||
//CGAL_Kernel_pred(Do_intersect_d,
|
||||
// do_intersect_d_object)
|
||||
CGAL_Kernel_pred(Less_lexicographically_d,
|
||||
less_lexicographically_d_object)
|
||||
CGAL_Kernel_pred(Less_or_equal_lexicographically_d,
|
||||
less_or_equal_lexicographically_d_object)
|
||||
CGAL_Kernel_pred(Linearly_independent_d,
|
||||
linearly_independent_d_object)
|
||||
CGAL_Kernel_pred(Linear_rank_d,
|
||||
linear_rank_d_object)
|
||||
CGAL_Kernel_pred(Orientation_d,
|
||||
orientation_d_object)
|
||||
CGAL_Kernel_pred(Coaffine_orientation_d,
|
||||
coaffine_orientation_d_object)
|
||||
CGAL_Kernel_pred(Side_of_bounded_sphere_d,
|
||||
side_of_bounded_sphere_d_object)
|
||||
CGAL_Kernel_pred(Side_of_oriented_sphere_d,
|
||||
side_of_oriented_sphere_d_object)
|
||||
CGAL_Kernel_pred(Side_of_oriented_subsphere_d,
|
||||
side_of_oriented_subsphere_d_object)
|
||||
CGAL_Kernel_pred(Oriented_side_d,
|
||||
oriented_side_d_object)
|
||||
|
||||
|
||||
CGAL_Kernel_cons(Linear_base_d,
|
||||
linear_base_d_object)
|
||||
CGAL_Kernel_cons(Center_of_sphere_d,
|
||||
center_of_sphere_d_object)
|
||||
CGAL_Kernel_cons(Intersection_d_,
|
||||
intersection_d_object)
|
||||
CGAL_Kernel_cons(Lift_to_paraboloid_d,
|
||||
lift_to_paraboloid_d_object)
|
||||
CGAL_Kernel_cons(Midpoint_d,
|
||||
midpoint_d_object)
|
||||
CGAL_Kernel_cons(Project_along_d_axis_d,
|
||||
project_along_d_axis_d_object)
|
||||
CGAL_Kernel_cons(Squared_distance_d,
|
||||
squared_distance_d_object)
|
||||
|
||||
#undef CGAL_Kernel_pred
|
||||
#undef CGAL_Kernel_cons
|
||||
#undef CGAL_Kernel_obj
|
||||
|
|
@ -118,5 +118,6 @@ public:
|
|||
} //namespace CGAL
|
||||
|
||||
#include <CGAL/Kernel_d/Linear_algebraCd_impl.h>
|
||||
#include <CGAL/Kernel_d/Interval_linear_algebra.h>
|
||||
|
||||
#endif // CGAL_LINEAR_ALGEBRACD_H
|
||||
|
|
|
|||
|
|
@ -0,0 +1,19 @@
|
|||
|
||||
#ifndef REFERENCED_ARGUMENT_H
|
||||
#define REFERENCED_ARGUMENT_H
|
||||
|
||||
template< typename T>
|
||||
struct Referenced_argument
|
||||
{
|
||||
Referenced_argument() : arg_() {}
|
||||
template<typename A>
|
||||
Referenced_argument(A a) : arg_(a) {}
|
||||
template<typename A, typename B>
|
||||
Referenced_argument(A a, B b) : arg_(a, b) {}
|
||||
T & arg() { return arg_; }
|
||||
|
||||
protected:
|
||||
T arg_;
|
||||
};
|
||||
|
||||
#endif // REFERENCED_ARGUMENT_H
|
||||
|
|
@ -80,6 +80,30 @@ This is the sign of the determinant
|
|||
typedef typename PointD::R R;
|
||||
typename R::Orientation_d or_; return or_(first,last); }
|
||||
|
||||
|
||||
template <class ForwardIterator>
|
||||
Orientation
|
||||
coaffine_orientation(ForwardIterator first, ForwardIterator last)
|
||||
/*{\Mfunc determines the orientation of the points in the tuple
|
||||
|A = tuple [first,last)| where $A$ consists of $k + 1$ points in $d$ - space.
|
||||
This is the sign of the determinant
|
||||
\[ \left\Lvert \begin{array}{cccc}
|
||||
1 & 1 & 1 & 1 \\
|
||||
A[0] & A[1] & \dots & A[k]
|
||||
\end{array} \right\Lvert \]
|
||||
where |A[a_i]| denotes the cartesian coordinate vector of
|
||||
the $i$-th point in $A$, after projection on an axis-aligned |k|-affine subspace.
|
||||
\precond |size [first,last) <= d| and |A[i].dimension() == d|
|
||||
$\forall 0 \leq i \leq d$
|
||||
and the value type of |ForwardIterator| is |Point_d<R>|.}*/
|
||||
{
|
||||
typedef typename std::iterator_traits<ForwardIterator>::value_type PointD;
|
||||
typedef typename PointD::R R;
|
||||
typename R::Coaffine_orientation_d or_;
|
||||
return or_(first,last);
|
||||
}
|
||||
|
||||
|
||||
template <class R, class ForwardIterator>
|
||||
Oriented_side side_of_oriented_sphere(
|
||||
ForwardIterator first, ForwardIterator last,
|
||||
|
|
|
|||
Loading…
Reference in New Issue