mirror of https://github.com/CGAL/cgal
686 lines
24 KiB
Plaintext
686 lines
24 KiB
Plaintext
\documentclass[a4paper,11pt,twoside]{article}
|
|
\usepackage{Lweb}
|
|
\input{defs}
|
|
|
|
\begin{document}
|
|
\title{Predicates and Constructions on Geometric Objects}
|
|
\author{M. Seel}
|
|
\maketitle
|
|
|
|
\section{Constructions on Points}
|
|
|
|
Lift a point $p$ to the paraboloid of revolution in dimension
|
|
$d+1$. The $d+1$st coordinate of $p'$ in $d+1$-space is the
|
|
sum of squares of the components of $p$.
|
|
<<Lift_to_paraboloidCd function object>>=
|
|
template <class R>
|
|
struct Lift_to_paraboloidCd {
|
|
typedef typename R::Point_d Point_d;
|
|
typedef typename R::FT FT;
|
|
typedef typename R::LA LA;
|
|
|
|
Point_d operator()(const Point_d& p) const
|
|
{
|
|
int d = p.dimension();
|
|
typename LA::Vector h(d+1);
|
|
FT sum = 0;
|
|
for (int i = 0; i<d; i++) {
|
|
h[i] = p.cartesian(i);
|
|
sum += h[i]*h[i];
|
|
}
|
|
h[d] = sum;
|
|
return Point_d(d+1,h.begin(),h.end());
|
|
}
|
|
};
|
|
|
|
@ We just ignore the $d$ component.
|
|
<<Project_along_d_axisCd function object>>=
|
|
template <class R>
|
|
struct Project_along_d_axisCd {
|
|
typedef typename R::Point_d Point_d;
|
|
typedef typename R::FT FT;
|
|
|
|
Point_d operator()(const Point_d& p) const
|
|
{ return Point_d(p.dimension()-1,
|
|
p.cartesian_begin(),p.cartesian_end()-1); }
|
|
};
|
|
|
|
@ Midpoint can be trivially obtained via vector arithmetic.
|
|
<<MidpointCd function object>>=
|
|
template <class R>
|
|
struct MidpointCd {
|
|
typedef typename R::Point_d Point_d;
|
|
Point_d operator()(const Point_d& p, const Point_d& q) const
|
|
{ return Point_d(p + (q-p)/2); }
|
|
};
|
|
|
|
@ The center of a sphere defined by a range of points is only defined
|
|
if the set of defining points are legal. If the defining points are
|
|
all equal the sphere is trivial. So assume otherwise. Then the center
|
|
$c$ is the unique point with equal distance to all the defining
|
|
points. A point $c$ has equal distance to point $p_0$ and $p_i$ if it
|
|
lies on the perpendicual bisector of $p_d$ and $p_i$, i.e., if it
|
|
satisfies the plane equation $(p_i - p_d)^T c = (p_i - p_d) (p_i +
|
|
p_d)/2$. Note that $p_i - p_d$ is the normal vector of the bisector
|
|
hyperplane and $(p_i + p_d)/2$ is the midpoint of $p_d$ and $p_i$. The
|
|
equation above translates into the equation
|
|
\[ \sum_{0 \le j < d} 2* (p_{ij} - p_{dj})c_j =
|
|
\sum_{0 \le j < d} (p_{ij} - p_{dj})(p_{ij} + p_{dj})
|
|
\]
|
|
for the homogeneous coordinates of the points and the center. We may
|
|
tentatively assume that $c_d = 1$, solve the corresponding linear
|
|
system, and then define the center.
|
|
<<Center_of_sphereCd function object>>=
|
|
template <class R>
|
|
struct Center_of_sphereCd {
|
|
typedef typename R::Point_d Point_d;
|
|
typedef typename R::FT FT;
|
|
typedef typename R::LA LA;
|
|
template <class Forward_iterator>
|
|
Point_d operator()(Forward_iterator start, Forward_iterator end) const
|
|
{ CGAL_assertion(start!=end);
|
|
int d = start->dimension();
|
|
typename LA::Matrix M(d);
|
|
typename LA::Vector b(d);
|
|
Point_d pd = *start++;
|
|
for (int i = 0; i < d; i++) {
|
|
// we set up the equation for p_i
|
|
Point_d pi = *start++;
|
|
b[i] = 0;
|
|
for (int j = 0; j < d; j++) {
|
|
M(i,j) = FT(2)*(pi.cartesian(j) - pd.cartesian(j));
|
|
b[i] += (pi.cartesian(j) - pd.cartesian(j)) *
|
|
(pi.cartesian(j) + pd.cartesian(j));
|
|
}
|
|
}
|
|
FT D;
|
|
typename LA::Vector x;
|
|
LA::linear_solver(M,b,x,D);
|
|
return Point_d(d,x.begin(),x.end());
|
|
}
|
|
|
|
}; // Center_of_sphereCd
|
|
|
|
|
|
@ The sqared distance operations just uses the inner product:
|
|
<<Squared_distanceCd function object>>=
|
|
template <class R>
|
|
struct Squared_distanceCd {
|
|
typedef typename R::Point_d Point_d;
|
|
typedef typename R::Vector_d Vector_d;
|
|
typedef typename R::FT FT;
|
|
FT operator()(const Point_d& p, const Point_d& q) const
|
|
{ Vector_d v = p-q; return v.squared_length(); }
|
|
};
|
|
|
|
@ Finally a predicate which allows to determine affine dependence of
|
|
three points together with a quotient determining the position of a
|
|
point with respect to the line. The operation returns true iff point
|
|
|p| lies on the line through |s| and |t|. The operation provides also
|
|
a quotient $\lambda$ such that $p = |s| + \lambda * |(t-s)|$ in this
|
|
case. \precond $s \neq t$. We just calculate the $\lambda$ in the
|
|
equation $p = s + \lambda * (t-s)$. The calculation idea is that
|
|
$\lambda = (p-s)/(t-s)$ (component wise) and one of the components
|
|
$t_i-s_i \neq 0$ as they are different.
|
|
<<Position_on_lineCd function object>>=
|
|
template <class R>
|
|
struct Position_on_lineCd {
|
|
typedef typename R::Point_d Point_d;
|
|
typedef typename R::LA LA;
|
|
typedef typename R::FT FT;
|
|
|
|
bool operator()(const Point_d& p, const Point_d& s, const Point_d& t,
|
|
FT& l) const
|
|
{ int d = p.dimension();
|
|
CGAL_assertion_msg((d==s.dimension())&&(d==t.dimension()&& d>0),
|
|
"position_along_line: argument dimensions disagree.");
|
|
CGAL_assertion_msg((s!=t),
|
|
"Position_on_line_d: line defining points are equal.");
|
|
FT lnum = (p.cartesian(0) - s.cartesian(0));
|
|
FT lden = (t.cartesian(0) - s.cartesian(0));
|
|
FT num(lnum), den(lden), lnum_i, lden_i;
|
|
for (int i = 1; i < d; i++) {
|
|
lnum_i = (p.cartesian(i) - s.cartesian(i));
|
|
lden_i = (t.cartesian(i) - s.cartesian(i));
|
|
if (lnum*lden_i != lnum_i*lden) return false;
|
|
if (lden_i != FT(0)) { den = lden_i; num = lnum_i; }
|
|
}
|
|
l = num/den; return true;
|
|
}
|
|
};
|
|
|
|
@ \section{Predicates on Points}
|
|
|
|
For an iterator range |[first,last)| we define |S = set [first,last)|
|
|
as the ordered tuple $(|S[0]|,|S[1]|, \ldots |S[d-1]|)$ where
|
|
$|S[i]| = |*| |++|^{(i)}|first|$ (the element obtained by $i$ times
|
|
forwarding the iterator by operator |++| and then dereferencing it to
|
|
get the value to which it points). We write |d = size [first,last)|.
|
|
This extends the syntax of random access iterators to input iterators.
|
|
If we index the tuple as above then we require that
|
|
$|++|^{(d+1)}|first == last|$.
|
|
|
|
In the following we require the Iterators to be globally
|
|
of value type |Point_d<R>|. Also if we are handed over an iterator
|
|
range |[first,last)|, then all points in |S = set [first,last)|
|
|
have the same dimension |dim|.
|
|
|
|
The barycentric coordinates of a point $p \in R^d$ in a affine space
|
|
of dimension $k$ spanned by the points $p_1, \ldots ,p_k$ are the
|
|
rational numbers $\lambda_i$ with $\sum_{i=0}^k \lambda_i p_i = p$ and
|
|
$\sum_{i=0}^k \lambda_i = 1$. Obviously the above conditions can be
|
|
written as a linear system $Mx=b$, where $M$ is the matrix obtained by
|
|
the $k$-tupel of points as column vectors and $b$ as $p$ interpreted
|
|
as a vector.
|
|
<<Barycentric_coordinatesCd function object>>=
|
|
template <class R>
|
|
struct Barycentric_coordinatesCd {
|
|
typedef typename R::Point_d Point_d;
|
|
typedef typename R::LA LA;
|
|
typedef typename R::FT FT;
|
|
|
|
template <class ForwardIterator, class OutputIterator>
|
|
OutputIterator operator()(ForwardIterator first, ForwardIterator last,
|
|
const Point_d& p, OutputIterator result)
|
|
{ TUPLE_DIM_CHECK(first,last,Barycentric_coordinates_d);
|
|
int n = std::distance(first,last);
|
|
int d = p.dimension();
|
|
typename R::Affine_rank_d affine_rank;
|
|
CGAL_assertion(affine_rank(first,last)==d);
|
|
std::vector< Point_d > V(first,last);
|
|
typename LA::Matrix M(d+1,V.size());
|
|
typename LA::Vector b(d+1), x;
|
|
for (int i=0; i<d; ++i) {
|
|
for (int j=0; j<V.size(); ++j)
|
|
M(i,j)=V[j].cartesian(i);
|
|
b[i] = p.cartesian(i);
|
|
}
|
|
for (int j=0; j<V.size(); ++j)
|
|
M(d,j) = 1;
|
|
b[d]=1;
|
|
FT D;
|
|
LA::linear_solver(M,b,x,D);
|
|
for (int i=0; i < x.dimension(); ++result, ++i) {
|
|
*result= x[i];
|
|
}
|
|
return result;
|
|
}
|
|
};
|
|
|
|
|
|
@ We determine the orientation of the points in the set |A = set
|
|
[first,last)| where $A$ consists of $d + 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[d]
|
|
\end{array} \right\Lvert \]
|
|
where |A[i]| denotes the cartesian coordinate vector of the $i$-th
|
|
point in $A$. \precond |size [first,last) == d+1| and
|
|
|A[i].dimension() == d| $\forall 0 \leq i \leq d$ and the value type
|
|
of |ForwardIterator| is |Point_d<R>|.
|
|
|
|
We set up this matrix and return its determinant. Actually, it is more
|
|
convenient to transpose it and to make the first row the last. This
|
|
changes the sign if the number of rows is even, i.e., if $d$ is odd.
|
|
<<OrientationCd function object>>=
|
|
template <class R>
|
|
struct OrientationCd {
|
|
typedef typename R::Point_d Point_d;
|
|
typedef typename R::LA LA;
|
|
|
|
template <class ForwardIterator>
|
|
Orientation operator()(ForwardIterator first, ForwardIterator last)
|
|
{ TUPLE_DIM_CHECK(first,last,Orientation_d);
|
|
int d = std::distance(first,last);
|
|
// range contains d points of dimension d-1
|
|
CGAL_assertion_msg(first->dimension() == d-1,
|
|
"Orientation_d: needs first->dimension() + 1 many points.");
|
|
typename LA::Matrix M(d); // quadratic
|
|
for (int i = 0; i < d; ++first,++i) {
|
|
for (int j = 0; j < d-1; ++j)
|
|
M(i,j) = first->cartesian(j);
|
|
M(i,d-1) = 1;
|
|
}
|
|
int row_correction = ( (d % 2 == 0) ? -1 : +1 );
|
|
// we invert the sign if the row number is even i.e. d is odd
|
|
return Orientation(row_correction * LA::sign_of_determinant(M));
|
|
}
|
|
};
|
|
|
|
@ We determine to which side of the sphere $S$ defined by the points
|
|
in |A = set [first,last)| the point $x$ belongs, where $A$ consists of
|
|
$d + 1$ points in $d$ - space. The positive side is determined by the
|
|
positive sign of the determinant
|
|
\[
|
|
\left\Lvert \begin{array}{ccccc}
|
|
1 & 1 & 1 & 1 & 1\\
|
|
|lift(A[0])| & |lift(A[1])| & \dots & |lift(A[d])| & |lift(x)|
|
|
\end{array} \right\Lvert
|
|
\]
|
|
where for a point $p$ with cartesian coordinates $p_i$ we use
|
|
|lift(p)| to denote the $d + 1$-dimensional point with cartesian
|
|
coordinate vector $(p_0, \ldots,p_{d-1},\sum_{0 \le i < d}p_i^2)$. If
|
|
the points in $A$ are positively oriented then the positive side is
|
|
the inside of the sphere and the negative side is the outside of the
|
|
sphere. \precond value type of |ForwardIterator| is |Point_d<R>|.
|
|
|
|
The lifting map |lift| maps a point |p| onto the paraboloid of
|
|
revolution and the position of the point |x| with respect to the
|
|
oriented sphere defined by the points in $A$ is the same as the
|
|
orientation of the lifted points. In order to evaluate the determinant
|
|
we multiply each column by the square of the homogenizing
|
|
coordinate. This turns any column into $(h_{d}^2,h_0
|
|
h_{d},\ldots,h_{d-1} h_{d},\sum_{0 \le i < d} h_i^2)$, where the
|
|
$h_i$'s denote homogeneous coordinates. Thus we set up this matrix and
|
|
return the sign of its determinant (it is actually simpler to set up
|
|
the transposed matrix and so that's what we do).
|
|
<<Side_of_oriented_sphereCd function object>>=
|
|
template <class R>
|
|
struct Side_of_oriented_sphereCd {
|
|
typedef typename R::Point_d Point_d;
|
|
typedef typename R::LA LA;
|
|
typedef typename R::FT FT;
|
|
|
|
template <class ForwardIterator>
|
|
Oriented_side operator()(ForwardIterator first, ForwardIterator last,
|
|
const Point_d& x)
|
|
{
|
|
TUPLE_DIM_CHECK(first,last,Side_of_oriented_sphere_d);
|
|
int d = std::distance(first,last); // |A| contains |d| points
|
|
CGAL_assertion_msg((d-1 == first->dimension()),
|
|
"Side_of_oriented_sphere_d: needs first->dimension()+1 many input points.");
|
|
typename LA::Matrix M(d + 1);
|
|
for (int i = 0; i < d; ++first, ++i) {
|
|
FT Sum = 0;
|
|
M(i,0) = 1;
|
|
for (int j = 0; j < d-1; j++) {
|
|
FT cj = first->cartesian(j);
|
|
M(i,j + 1) = cj; Sum += cj*cj;
|
|
}
|
|
M(i,d) = Sum;
|
|
}
|
|
FT Sum = 0;
|
|
M(d,0) = 1;
|
|
for (int j = 0; j < d-1; j++) {
|
|
FT hj = x.cartesian(j);
|
|
M(d,j + 1) = hj; Sum += hj*hj;
|
|
}
|
|
M(d,d) = Sum;
|
|
return - LA::sign_of_determinant(M);
|
|
}
|
|
};
|
|
|
|
@ We determine whether the point |p| lies |ON_BOUNDED_SIDE|,
|
|
|ON_BOUNDARY|, or |ON_UNBOUNDED_SIDE| of the sphere defined by
|
|
the points in |A = set [first,last)| where $A$ consists of $d+1$
|
|
points in $d$-space.
|
|
\precond value type of |ForwardIterator| is |Point_d<R>| and
|
|
$|orientation(first,last)| \neq |ZERO|$.
|
|
<<Side_of_bounded_sphereCd function object>>=
|
|
template <class R>
|
|
struct Side_of_bounded_sphereCd {
|
|
typedef typename R::Point_d Point_d;
|
|
|
|
template <class ForwardIterator>
|
|
Bounded_side operator()(ForwardIterator first, ForwardIterator last,
|
|
const Point_d& p)
|
|
{
|
|
TUPLE_DIM_CHECK(first,last,region_of_sphere);
|
|
typename R::Orientation_d _orientation;
|
|
Orientation o = _orientation(first,last);
|
|
CGAL_assertion_msg((o != 0), "Side_of_bounded_sphere_d: \
|
|
A must be full dimensional.");
|
|
typename R::Side_of_oriented_sphere_d _side_of_oriented_sphere;
|
|
Oriented_side oside = _side_of_oriented_sphere(first,last,p);
|
|
if (o == POSITIVE) {
|
|
switch (oside) {
|
|
case ON_POSITIVE_SIDE : return ON_BOUNDED_SIDE;
|
|
case ON_ORIENTED_BOUNDARY: return ON_BOUNDARY;
|
|
case ON_NEGATIVE_SIDE : return ON_UNBOUNDED_SIDE;
|
|
}
|
|
} else {
|
|
switch (oside) {
|
|
case ON_POSITIVE_SIDE : return ON_UNBOUNDED_SIDE;
|
|
case ON_ORIENTED_BOUNDARY: return ON_BOUNDARY;
|
|
case ON_NEGATIVE_SIDE : return ON_BOUNDED_SIDE;
|
|
}
|
|
}
|
|
return ON_BOUNDARY; // never reached
|
|
}
|
|
};
|
|
|
|
|
|
@ We determine whether |p| is contained in the simplex spanned by the
|
|
points in |A = set[first,last)|. |A| may consists of up to $d + 1$
|
|
points. \precond value type of |ForwardIterator| is |Point_d<R>| and
|
|
the points in |A| are affinely independent.
|
|
|
|
A point |p| is contained in the convex hull of a set |A| of points if
|
|
|p| is a convex combination of the points in |A|, i.e., if the system
|
|
$\sum \lambda_i A_i = p$ has a solution with $\sum \lambda_i = 1$ and
|
|
$\lambda_i \ge 0$. If the points in $A$ are linearly independent then
|
|
the solution is unique (if it exists at all). We therefore proceed as
|
|
above and then check whether the solution vector is non-negative.
|
|
<<Contained_in_simplexCd function object>>=
|
|
template <class R>
|
|
struct Contained_in_simplexCd {
|
|
typedef typename R::Point_d Point_d;
|
|
typedef typename R::LA LA;
|
|
typedef typename R::FT FT;
|
|
|
|
template <class ForwardIterator>
|
|
bool operator()(ForwardIterator first, ForwardIterator last,
|
|
const Point_d& p)
|
|
{
|
|
TUPLE_DIM_CHECK(first,last,Contained_in_simplex_d);
|
|
int k = std::distance(first,last); // |A| contains |k| points
|
|
int d = first->dimension();
|
|
CGAL_assertion_code(
|
|
typename R::Affinely_independent_d check_independence; )
|
|
CGAL_assertion_msg(check_independence(first,last),
|
|
"Contained_in_simplex_d: A not affinely independent.");
|
|
CGAL_assertion(d==p.dimension());
|
|
|
|
typename LA::Matrix M(d + 1,k);
|
|
typename LA::Vector b(d +1);
|
|
for (int j = 0; j < k; ++first, ++j) {
|
|
for (int i = 0; i < d; ++i)
|
|
M(i,j) = first->cartesian(i);
|
|
M(d,j) = 1;
|
|
}
|
|
for (int i = 0; i < d; ++i)
|
|
b[i] = p.cartesian(i);
|
|
b[d] = 1;
|
|
|
|
FT D;
|
|
typename LA::Vector lambda;
|
|
if ( LA::linear_solver(M,b,lambda,D) ) {
|
|
for (int j = 0; j < k; j++) {
|
|
if (lambda[j] < FT(0)) return false;
|
|
}
|
|
return true;
|
|
}
|
|
return false;
|
|
}
|
|
};
|
|
|
|
@ We determine whether $p$ is contained in the affine hull of the
|
|
points in |A = set[first,last)|. \precond value type of
|
|
|ForwardIterator| is |Point_d<R>|.
|
|
|
|
A point |p| is contained in the affine hull of a set |A| of points if
|
|
|p| is an affine combination of the points in |A|, i.e., if the system
|
|
$\sum \lambda_i A_i = p$ has a solution with $\sum \lambda_i = 1$. Set
|
|
$\lambda_i = A_{i,d} \beta_i /p_d$ with $A_{i,d}$ being the
|
|
homogenizing coordinate of $A_i$ and $p_d$ being the homogenizing
|
|
coordinate of $p$. The $i$-th column of the system for the
|
|
$\beta_i$'s is simply the homogeneous vector of $A_i$ and the right
|
|
hand side is simply the homogeneous vector for $p$. Thus we proceed as
|
|
above but let $i$ run up to $d$.
|
|
<<Contained_in_affine_hullCd function object>>=
|
|
template <class R>
|
|
struct Contained_in_affine_hullCd {
|
|
typedef typename R::Point_d Point_d;
|
|
typedef typename R::LA LA;
|
|
|
|
template <class ForwardIterator>
|
|
bool operator()(ForwardIterator first, ForwardIterator last,
|
|
const Point_d& p)
|
|
{
|
|
TUPLE_DIM_CHECK(first,last,Contained_in_affine_hull_d);
|
|
int k = std::distance(first,last); // |A| contains |k| points
|
|
int d = first->dimension();
|
|
typename LA::Matrix M(d + 1,k);
|
|
typename LA::Vector b(d + 1);
|
|
for (int j = 0; j < k; ++first, ++j) {
|
|
for (int i = 0; i < d; ++i)
|
|
M(i,j) = first->cartesian(i);
|
|
M(d,j) = 1;
|
|
}
|
|
for (int i = 0; i < d; ++i)
|
|
b[i] = p.cartesian(i);
|
|
b[d] = 1;
|
|
return LA::is_solvable(M,b);
|
|
}
|
|
};
|
|
|
|
|
|
@ We compute the affine rank of the points in |A = set [first,last)|.
|
|
\precond value type of |ForwardIterator| is |Point_d<R>|.
|
|
|
|
The affine rank of the $k+1$ points $p_0$, $p_1$,\ldots, $p_k$ is the
|
|
linear rank of the $k$ vectors $p_1 - p_0$, \ldots, $p_k - p_0$.
|
|
<<Affine_rankCd function object>>=
|
|
template <class R>
|
|
struct Affine_rankCd {
|
|
typedef typename R::Point_d Point_d;
|
|
typedef typename R::Vector_d Vector_d;
|
|
typedef typename R::LA LA;
|
|
|
|
template <class ForwardIterator>
|
|
int operator()(ForwardIterator first, ForwardIterator last)
|
|
{
|
|
TUPLE_DIM_CHECK(first,last,Affine_rank_d);
|
|
int k = std::distance(first,last); // |A| contains |k| points
|
|
if (k == 0) return -1;
|
|
if (k == 1) return 0;
|
|
int d = first->dimension();
|
|
typename LA::Matrix M(d,--k);
|
|
Point_d p0 = *first; ++first; // first points to second
|
|
for (int j = 0; j < k; ++first, ++j) {
|
|
Vector_d v = *first - p0;
|
|
for (int i = 0; i < d; i++)
|
|
M(i,j) = v.cartesian(i);
|
|
}
|
|
return LA::rank(M);
|
|
}
|
|
};
|
|
|
|
@ We decide whether the points in |A = set [first,last)| are affinely
|
|
independent. \precond value type of |ForwardIterator| is |Point_d<R>|.
|
|
A set of points is affinely independent if their affine rank is equal
|
|
to the number of points minus 1.
|
|
<<Affinely_independentCd function object>>=
|
|
template <class R>
|
|
struct Affinely_independentCd {
|
|
typedef typename R::Point_d Point_d;
|
|
typedef typename R::LA LA;
|
|
|
|
template <class ForwardIterator>
|
|
bool operator()(ForwardIterator first, ForwardIterator last)
|
|
{ typename R::Affine_rank_d rank;
|
|
int n = std::distance(first,last);
|
|
return rank(first,last) == n-1;
|
|
}
|
|
};
|
|
|
|
|
|
@ We compare the Cartesian coordiantes of points |p1| and |p2|
|
|
lexicographically.
|
|
<<Compare_lexicographicallyCd function object>>=
|
|
template <class R>
|
|
struct Compare_lexicographicallyCd {
|
|
typedef typename R::Point_d Point_d;
|
|
typedef typename R::Point_d PointD; //MSVC hack
|
|
Comparison_result operator()(const Point_d& p1, const Point_d& p2)
|
|
{ return PointD::cmp(p1,p2); }
|
|
};
|
|
|
|
@ \section{Predicates on Vectors}
|
|
|
|
We determine whether $x$ is contained in the linear hull of the
|
|
vectors in |A = set[first,last)|. \precond value type of
|
|
|ForwardIterator| is |Vector_d<R>|.
|
|
|
|
A vector |x| is contained in the linear hull of a set |A| of vectors
|
|
if |x| is a linear combination of the vectors in |A|, i.e., if the
|
|
system $\sum \lambda_i A_i = x$ has a solution.
|
|
<<Contained_in_linear_hullCd function object>>=
|
|
template <class R>
|
|
struct Contained_in_linear_hullCd {
|
|
typedef typename R::LA LA;
|
|
typedef typename R::FT FT;
|
|
typedef typename R::Vector_d Vector_d;
|
|
|
|
template<class ForwardIterator>
|
|
bool operator()(
|
|
ForwardIterator first, ForwardIterator last, const Vector_d& x)
|
|
{ TUPLE_DIM_CHECK(first,last,Contained_in_linear_hull_d);
|
|
int k = std::distance(first,last); // |A| contains |k| vectors
|
|
int d = first->dimension();
|
|
typename LA::Matrix M(d,k);
|
|
typename LA::Vector b(d);
|
|
for (int i = 0; i < d; i++) {
|
|
b[i] = x.cartesian(i);
|
|
for (int j = 0; j < k; j++)
|
|
M(i,j) = (first+j)->cartesian(i);
|
|
}
|
|
return LA::is_solvable(M,b);
|
|
}
|
|
};
|
|
|
|
@ We compute the linear rank of the vectors in |A = set [first,last)|.
|
|
\precond value type of |ForwardIterator| is |Vector_d<R>|.
|
|
|
|
We set up a matrix having the cartesian coordinates of the vectors in
|
|
A as its columns. The linear rank is the rank of this matrix.
|
|
<<Linear_rankCd function object>>=
|
|
template <class R>
|
|
struct Linear_rankCd {
|
|
typedef typename R::LA LA;
|
|
template <class ForwardIterator>
|
|
int operator()(ForwardIterator first, ForwardIterator last)
|
|
{ TUPLE_DIM_CHECK(first,last,linear_rank);
|
|
int k = std::distance(first,last); // k vectors
|
|
int d = first->dimension();
|
|
typename LA::Matrix M(d,k);
|
|
for (int i = 0; i < d ; i++)
|
|
for (int j = 0; j < k; j++)
|
|
M(i,j) = (first + j)->cartesian(i);
|
|
return LA::rank(M);
|
|
}
|
|
};
|
|
|
|
@ We decide whether the vectors in $A$ are linearly independent.
|
|
\precond value type of |ForwardIterator| is |Vector_d<R>|.
|
|
|
|
A set of vectors is linearly independent if their linear rank is
|
|
equal to the number of vectors in the set.
|
|
<<Linearly_independentCd function object>>=
|
|
template <class R>
|
|
struct Linearly_independentCd {
|
|
typedef typename R::LA LA;
|
|
template <class ForwardIterator>
|
|
bool operator()(ForwardIterator first, ForwardIterator last)
|
|
{ typename R::Linear_rank_d rank;
|
|
return rank(first,last) == std::distance(first,last);
|
|
}
|
|
};
|
|
|
|
@ We compute a basis of the linear space spanned by the vectors in
|
|
|set [first,last)| and return it via an iterator range starting in
|
|
|result|. The returned iterator marks the end of the output. \precond
|
|
value type of |ForwardIterator| and |OutputIterator| is
|
|
|Vector_d<R>|. To find a basis within |A| we use the corresponding
|
|
function |independent_columns| of our linear algebra package.
|
|
<<Linear_baseCd function object>>=
|
|
template <class R>
|
|
struct Linear_baseCd {
|
|
typedef typename R::LA LA;
|
|
typedef typename R::FT FT;
|
|
typedef typename R::Vector_d Vector_d;
|
|
template <class ForwardIterator, class OutputIterator>
|
|
OutputIterator operator()(ForwardIterator first, ForwardIterator last,
|
|
OutputIterator result)
|
|
{ TUPLE_DIM_CHECK(first,last,linear_base);
|
|
int k = std::distance(first,last); // k vectors
|
|
int d = first->dimension();
|
|
typename LA::Matrix M(d,k);
|
|
for (int j = 0; j < k; ++first, ++j)
|
|
for (int i = 0; i < d; i++)
|
|
M(i,j) = first->cartesian(i);
|
|
|
|
std::vector<int> indcols;
|
|
int r = LA::independent_columns(M,indcols);
|
|
|
|
for (int l=0; l < r; l++) {
|
|
typename LA::Vector v = M.column(indcols[l]);
|
|
*result++ = Vector_d(d,v.begin(),v.end());
|
|
}
|
|
return result;
|
|
}
|
|
};
|
|
|
|
|
|
<<function_objectsCd.h>>=
|
|
// ======================================================================
|
|
//
|
|
// Copyright (c) 2000,2001 The CGAL Consortium
|
|
//
|
|
// This software and related documentation is part of an INTERNAL release
|
|
// of the Computational Geometry Algorithms Library (CGAL). It is not
|
|
// intended for general use.
|
|
//
|
|
// ----------------------------------------------------------------------
|
|
//
|
|
// release : $CGAL_Revision$
|
|
// release_date : $CGAL_Date$
|
|
//
|
|
// file : include/CGAL/Kernel_d/function_objectsCd.h
|
|
// package : Kernel_d
|
|
// maintainer : Michael Seel <Michael.Seel@mpi-sb.mpg.de>
|
|
// revision : $Id$
|
|
// revision_date : $Date$
|
|
// author(s) : Michael Seel
|
|
// coordinator : MPI Saarbruecken (Susan.Hert@mpi-sb.mpg.de)
|
|
//
|
|
// ======================================================================
|
|
//---------------------------------------------------------------------
|
|
// file generated by notangle from noweb/function_objectsCd.lw
|
|
// please debug or modify noweb file
|
|
// coding: K. Mehlhorn, M. Seel
|
|
//---------------------------------------------------------------------
|
|
#ifndef CGAL_FUNCTION_OBJECTSCD_H
|
|
#define CGAL_FUNCTION_OBJECTSCD_H
|
|
|
|
#ifndef NOCGALINCL
|
|
#include <CGAL/basic.h>
|
|
#include <CGAL/enum.h>
|
|
#endif
|
|
|
|
#undef TRACE
|
|
#undef TRACEN
|
|
#undef TRACEV
|
|
#define TRACE(t) std::cerr << t
|
|
#define TRACEN(t) std::cerr << t << std::endl
|
|
#define TRACEV(t) std::cerr << #t << " = " << (t) << std::endl
|
|
|
|
CGAL_BEGIN_NAMESPACE
|
|
|
|
<<Lift_to_paraboloidCd function object>>
|
|
<<Project_along_d_axisCd function object>>
|
|
<<MidpointCd function object>>
|
|
<<Center_of_sphereCd function object>>
|
|
<<Squared_distanceCd function object>>
|
|
<<Position_on_lineCd function object>>
|
|
<<Barycentric_coordinatesCd function object>>
|
|
<<OrientationCd function object>>
|
|
<<Side_of_oriented_sphereCd function object>>
|
|
<<Side_of_bounded_sphereCd function object>>
|
|
<<Contained_in_simplexCd function object>>
|
|
<<Contained_in_affine_hullCd function object>>
|
|
<<Affine_rankCd function object>>
|
|
<<Affinely_independentCd function object>>
|
|
<<Compare_lexicographicallyCd function object>>
|
|
<<Contained_in_linear_hullCd function object>>
|
|
<<Linear_rankCd function object>>
|
|
<<Linearly_independentCd function object>>
|
|
<<Linear_baseCd function object>>
|
|
|
|
CGAL_END_NAMESPACE
|
|
#endif //CGAL_FUNCTION_OBJECTSCD_H
|
|
|
|
@
|
|
\end{document}
|