mirror of https://github.com/CGAL/cgal
added missing files
This commit is contained in:
Michael Seel
parent
19b10d368c
commit
0f10a4e6af
|
|
@ -0,0 +1 @@
|
|||
noweb
|
||||
|
|
@ -0,0 +1,662 @@
|
|||
\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_circleCd function object>>=
|
||||
template <class R>
|
||||
struct Center_of_circleCd {
|
||||
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_circleCd
|
||||
|
||||
|
||||
@ 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 != 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 (register int i=0; i<d; ++i) {
|
||||
for (register int j=0; j<V.size(); ++j)
|
||||
M(i,j)=V[j].cartesian(i);
|
||||
b[i] = p.cartesian(i);
|
||||
}
|
||||
for (register 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 (register int i = 0; i < d; ++first, ++i) {
|
||||
FT Sum = 0;
|
||||
M(i,0) = 1;
|
||||
for (register 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 (register 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 Oriented_side( - 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 or = _orientation(first,last);
|
||||
CGAL_assertion_msg((or != 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 (or == 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();
|
||||
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 (register int j = 0; j < k; ++first, ++j) {
|
||||
for (register int i = 0; i < d; ++i)
|
||||
M(i,j) = first->cartesian(i);
|
||||
M(d,j) = 1;
|
||||
}
|
||||
for (register 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] < 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 (register int j = 0; j < k; ++first, ++j) {
|
||||
for (register int i = 0; i < d; ++i)
|
||||
M(i,j) = first->cartesian(i);
|
||||
M(d,j) = 1;
|
||||
}
|
||||
for (register 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;
|
||||
Comparison_result operator()(const Point_d& p1, const Point_d& p2)
|
||||
{ return Point_d::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();
|
||||
FT denom;
|
||||
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>>=
|
||||
//---------------------------------------------------------------------
|
||||
// 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_circleCd 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}
|
||||
Loading…
Reference in New Issue