songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity
cgal/Packages/Kernel_d/noweb/HyperplaneHd.lw

504 lines
18 KiB
Plaintext
Raw Blame History

\documentclass[a4paper,11pt,twoside]{article}
\usepackage{Lweb}
\begin{document}
\title{Hyperplanes with Rational Coordinates in d-Space\\
(class HyperplaneHd)}
\author{M. Seel}
\date{\today}
\maketitle
\tableofcontents
\newpage
@ \section{The Manual Page of class HyperplaneHd}
\input HyperplaneHd.man
@ \section{The Implementation of class HyperplaneHd}
The type HyperplaneHd is an item class with representation class georep. It
shares this representation class with points, vectors, and
directions. We derive HyperplaneHd from handle\_base and derive georep from
handle\_rep. This gives us reference counting for free. We give all
implementations which are trivial directly in the header file and
postpone all others to the next section. Aside from this the header
file is in one-to-one correspondence to the manual page.
<<HyperplaneHd.h>>=
//---------------------------------------------------------------------
// file generated by notangle from HyperplaneHd.lw
// please debug or modify noweb file
// coding: K. Mehlhorn, M. Seel
//---------------------------------------------------------------------
#ifndef CGAL_HYPERPLANEHD_H
#define CGAL_HYPERPLANEHD_H
#ifndef NOCGALINCL
#include <CGAL/basic.h>
#include <CGAL/Quotient.h>
#endif
#include <CGAL/Kernel_d/PointHd.h>
#include <CGAL/Kernel_d/VectorHd.h>
#include <CGAL/Kernel_d/Aff_transformationHd.h>
CGAL_BEGIN_NAMESPACE
<<prototyping>>
<<defining HyperplaneHd>>
CGAL_END_NAMESPACE
#endif // CGAL_HYPERPLANEHD_H
<<HyperplaneHd.C>>=
//---------------------------------------------------------------------
// file generated by notangle from HyperplaneHd.lw
// please debug or modify noweb file
// coding: K. Mehlhorn, M. Seel
//---------------------------------------------------------------------
#ifndef CGAL_HYPERPLANEHD_C
#define CGAL_HYPERPLANEHD_C
CGAL_BEGIN_NAMESPACE
<<implementing HyperplaneHd>>
CGAL_END_NAMESPACE
#endif // CGAL_HYPERPLANEHD_C
@ \subsection{The Datatype}
And now for the class definition. We interleave the prototyping and
the implementation.
<<defining HyperplaneHd>>=
/*{\Manpage{Hyperplane_d}{R}{Hyperplanes in d-space}{h}}*/
/*{\Msubst
Hd<RT,LA>#_d<R>
HyperplaneHd#Hyperplane_d
Quotient<RT>#FT
}*/
template <class _RT, class _LA>
class HyperplaneHd : public Handle_for< Tuple_d<_RT,_LA> > {
typedef Tuple_d<_RT,_LA> Tuple;
typedef Handle_for<Tuple> Base;
typedef HyperplaneHd<_RT,_LA> Self;
/*{\Mdefinition An instance of data type |HyperplaneHd| is an
oriented hyperplane in $d$ - dimensional space. A hyperplane $h$ is
represented by coefficients $(c_0,c_1,\ldots,c_d)$ of type |RT|. At
least one of $c_0$ to $c_{ d - 1 }$ must be non-zero. The plane
equation is $\sum_{ 0 \le i < d } c_i x_i + c_d = 0$, where $x_0$ to
$x_{d-1}$ are Cartesian point coordinates.
For a particular $x$ the sign of $\sum_{ 0 \le i < d } c_i x_i +
c_d$ determines the position of a point $x$ with respect to the
hyperplane (on the hyperplane, on the negative side, or on the
positive side).
There are two equality predicates for hyperplanes. The (weak)
equality predicate (|weak_equality|) declares two hyperplanes equal if
they consist of the same set of points, the strong equality predicate
(|operator==|) requires in addition that the negative halfspaces
agree. In other words, two hyperplanes are strongly equal if their
coefficient vectors are positive multiples of each other and they are
(weakly) equal if their coefficient vectors are multiples of each
other.}*/
const typename _LA::Vector& vector_rep() const { return ptr->v; }
_RT& entry(int i) const { return ptr->v[i]; }
void invert_rep() { ptr->invert(); }
public:
/*{\Mtypes 4}*/
typedef _RT RT;
/*{\Mtypemember the ring type.}*/
typedef Quotient<_RT> FT;
/*{\Mtypemember the field type.}*/
typedef _LA LA;
/*{\Mtypemember the linear algebra layer.}*/
typedef typename Tuple::const_iterator Coefficient_const_iterator;
/*{\Mtypemember a read-only iterator for the coefficients.}*/
/*{\Mcreation h 4}*/
/*{\Moptions nextwarning=no}*/
HyperplaneHd(int d = 0) : Base( Tuple(d+1) ) {}
/*{\Mcreate introduces a variable |\Mvar| of type |\Mname|
initialized to some hyperplane in $d$ - dimensional space. }*/
template <class InputIterator>
HyperplaneHd(int d, InputIterator first, InputIterator last, RT D)
: Base( Tuple(d+1,first,last,D) ) {}
/*{\Mcreate introduces a variable |\Mvar| of type |\Mname|
initialized to the hyperplane with coefficients |set [first,last)| and
|D|. \precond |size [first,last) == d| and the value type of
InputIterator is |RT|.}*/
template <class InputIterator>
HyperplaneHd(int d, InputIterator first, InputIterator last)
: Base( Tuple(d+1,first,last) ) {}
/*{\Mcreate introduces a variable |\Mvar| of type |\Mname|
initialized to the hyperplane with coefficients |set [first,last)|.
\precond |size [first,last) == d+1| and the value type of
InputIterator is |RT|.}*/
template <class ForwardIterator>
HyperplaneHd(ForwardIterator first, ForwardIterator last,
const PointHd<RT,LA>& o,
Oriented_side side = Oriented_side(0))
/*{\Mcreate constructs some hyperplane that passes through the points
in |set [first,last)|. If |side| is |ON_POSITIVE_SIDE| or
|ON_NEGATIVE_SIDE| then |o| is on that side of the constructed
hyperplane. \precond A hyperplane with the stated properties must
exist. The value type of |ForwardIterator| is |PointHd<RT,LA>|. }*/
@ We want to construct a hyperplane that passes through a set |P = set
[first,last)| of points in $d$-dimensional space and has a specified
point $o$ on a specified side. We simply have to find a vector $x$
such that $P^T \cdot x = 0$ for every point in $P$. This amounts to
solving a homogeneous linear system. If the system has only a trivial
solution the task at hand is unsolvable and we report an error. So
assume that the system has a non-trivial solution. Let vectors $s_1,
\ldots, s_k$ span the solution space. if |side == ZERO| we may take
any $s_j$ as the normal vector of our hyperplane. if $|side| \neq 0$
and the task at hand is solvable there must be a $j$ such that $o^T
\cdot s_j \neq 0$. We take $s_j$ as the normal vector of our
hyperplane and use |o| to normalize the hyperplane equation.
<<defining HyperplaneHd>>=
: Base( Tuple(o.dimension()+1) ) {
TUPLE_DIM_CHECK(first,last,hyperplane::construction);
CGAL_assertion_msg((first->dimension()==o.dimension()),
"hyperplane::construction: dimensions disagree.");
int d = first->dimension(); // we are in $d$ - dimensional space
int m = std::distance(first,last); // |P| has $m$ points
typename LA::Matrix A(m,d + 1);
for (int i = 0; i < m; i++) { /* define $i$-th equation */
for (int j = 0; j <= d; j++)
A(i,j) = first->homogeneous(j); // $j$ - th coord of $i$-th point
++first;
}
typename LA::Matrix spanning_vecs; // columns span solution
int dim = LA::homogeneous_linear_solver(A,spanning_vecs);
if (dim == 0)
CGAL_assertion_msg(0,"HyperplaneHd::constructor: \
set P is full dimensional.");
if (side == ON_ORIENTED_BOUNDARY) {
ptr->v = spanning_vecs.column(0);
return;
}
RT sum = 0;
int j;
for (j = 0; j < dim; j++) {
for (int i = 0; i <= d; i++)
sum += spanning_vecs(i,j)*o.homogeneous(i);
if (sum != 0) break;
}
if (j == dim)
CGAL_assertion_msg(0,"HyperplaneHd::constructor: \
cannot use o to determine side.");
ptr->v = spanning_vecs.column(j);
if ( CGAL_NTS sign(sum) > 0 && side == ON_NEGATIVE_SIDE ||
CGAL_NTS sign(sum) < 0 && side == ON_POSITIVE_SIDE)
invert_rep();
}
HyperplaneHd(const PointHd<RT,LA>& p, const DirectionHd<RT,LA>& dir)
/*{\Mcreate constructs the hyperplane with normal direction |dir|
that passes through $p$. The direction |dir| points into the positive
side. \precond |dir| is not the trivial direction.}*/
@ Given a point |p| and a direction |dir| we want to construct a hyperplane
with normal direction |dir| and passing through |p|. We set the coefficient
vector $x = (|dir|_0, \ldots,|dir|_{d-1},D)$ for some unknown $D$ and then
use |p| to determine $D$ such that $p^T \cdot x = 0$.
Note that $D$ will be rational in general.
<<defining HyperplaneHd>>=
: Base( Tuple(p.dimension()+1) ) {
int d = p.dimension();
CGAL_assertion_msg((dir.dimension() == d), "HyperplaneHd::constructor: \
parameter dimensions disagree.");
CGAL_assertion_msg((dir.dimension() == d), "HyperplaneHd::constructor: \
parameter dimensions disagree.");
RT sum = 0;
for (int i = 0; i < d; i++) {
sum += dir.delta(i)*p.homogeneous(i);
entry(i) = dir.delta(i)*p.homogeneous(d);
}
entry(d) = -sum;
}
HyperplaneHd(const RT& a, const RT& b, const RT& c) :
Base( Tuple(a,b,c) ) {}
/*{\Mcreate introduces a variable |\Mvar| of type |\Mname| in
$2$-dimensional space with equation $ax+by+c=0$. }*/
HyperplaneHd(int a, int b, int c) :
Base( Tuple(RT(a),RT(b),RT(c)) ) {}
HyperplaneHd(const RT& a, const RT& b, const RT& c, const RT& d) :
Base( Tuple(a,b,c,d) ) {}
/*{\Mcreate introduces a variable |\Mvar| of type |\Mname| in
$3$-dimensional space with equation $ax+by+cz+d=0$. }*/
HyperplaneHd(int a, int b, int c, int d) :
Base( Tuple(RT(a),RT(b),RT(c),RT(d)) ) {}
HyperplaneHd(const HyperplaneHd<RT,LA>& h) : Base(h) {}
~HyperplaneHd() {}
/*{\Moperations 4 2}*/
int dimension() const { return ptr->size()-1; }
/*{\Mop returns the dimension of |\Mvar|. }*/
RT operator[](int i) const
/*{\Marrop returns the $i$-th coefficient of |\Mvar|.
\precond $0 \leq i \leq d$.}*/
{ CGAL_assertion_msg((0<=i && i<=(dimension())), "HyperplaneHd::op[]:\
index out of range."); return entry(i); }
RT coefficient(int i) const { return entry(i); }
/*{\Mop returns the $i$-th coefficient of |\Mvar|.
\precond $0 \leq i \leq d$.}*/
const typename LA::Vector& coefficient_vector() const
/*{\Xop returns the coefficient vector $(c_0,\ldots,c_d)$ of |\Mvar|. }*/
{ return vector_rep(); }
Coefficient_const_iterator coefficients_begin() const
/*{\Mop returns an iterator pointing to the first coefficient.}*/
{ return ptr->begin(); }
Coefficient_const_iterator coefficients_end() const
/*{\Mop returns an iterator pointing beyond the last coefficient.}*/
{ return ptr->end(); }
VectorHd<RT,LA> orthogonal_vector() const;
/*{\Mop returns the orthogonal vector of |\Mvar|. It points from the
negative halfspace into the positive halfspace and its
homogeneous coordinates are $(c_0, \ldots, c_{d - 1},1)$. }*/
@ Any multiple of $(c_0,\ldots,c_{d-1})$ is a normal vector. We want the
vector to point from the negative to the positive halfspace. Recall
that our hyperplane has the equation $c_d + \sum (c_i z_i) = 0$, where
the $z_i$ are Cartesian point coordinates. The point $z = -c_d\cdot
c/{\parallel c \parallel}$ is on the hyperplane, the point $z_n = (-1
-c_d)\cdot c/{\parallel c \parallel}$ is in the negative halfspace and
the point $z_p = (1 - c_d)\cdot c/{\parallel c \parallel}$ is in the
positive halfspace. Thus any positive multiple of $c$ is the desired
orthogonal vector. We take $(c_0,\ldots,c_{d-1},1)$.
<<implementing HyperplaneHd>>=
template <class RT, class LA>
VectorHd<RT,LA> HyperplaneHd<RT,LA>::
orthogonal_vector() const
{ VectorHd<RT,LA> res(*this);
res.copy_on_write();
res.entry(dimension()) = 1;
return res;
}
<<defining HyperplaneHd>>=
DirectionHd<RT,LA> orthogonal_direction() const
/*{\Mop returns the orthogonal direction of |\Mvar|. It points from the
negative halfspace into the positive halfspace. }*/
{ return orthogonal_vector().direction(); }
RT value_at(const PointHd<RT,LA>& p) const
/*{\Xop returns the value of |\Mvar| at the point |p|, i.e.,
$\sum_{ 0 \le i \le d } h_i p_i$.\\
Warning: this value depends on the particular representation
of |\Mvar| and |p|. }*/
{ CGAL_assertion_msg((dimension()==p.dimension()),"HyperplaneHd::value_at:\
dimensions disagree.");
return vector_rep()*p.vector_rep();
}
Oriented_side oriented_side(const PointHd<RT,LA>& p) const
/*{\Mop returns the side of the hyperplane |\Mvar| containing $p$. }*/
/*{\Mtext \setopdims{2cm}{2cm}}*/
{
CGAL_assertion_msg((dimension()==p.dimension()),
"HyperplaneHd::oriented_side: dimensions do not agree.");
return Oriented_side(CGAL_NTS sign(value_at(p)));
}
bool has_on(const PointHd<RT,LA>& p) const
/*{\Mop returns true iff point |p| lies on the hyperplane |\Mvar|. }*/
{ return (oriented_side(p) == ON_ORIENTED_BOUNDARY); }
bool has_on_boundary(const PointHd<RT,LA>& p) const
/*{\Mop returns true iff point |p| lies on the boundary of
hyperplane |\Mvar|. }*/
{ return (oriented_side(p) == ON_ORIENTED_BOUNDARY); }
bool has_on_positive_side(const PointHd<RT,LA>& p) const
/*{\Mop returns true iff point |p| lies on the positive side of
hyperplane |\Mvar|. }*/
{ return (oriented_side(p) == ON_POSITIVE_SIDE); }
bool has_on_negative_side(const PointHd<RT,LA>& p) const
/*{\Mop returns true iff point |p| lies on the negative side of
hyperplane |\Mvar|. }*/
{ return (oriented_side(p) == ON_NEGATIVE_SIDE); }
/*{\Mtext \restoreopdims }*/
HyperplaneHd<RT,LA> transform(const Aff_transformationHd<RT,LA>& t) const
/*{\Mop returns $t(h)$.}*/
{ typename LA::Vector res =
-(LA::transpose(t.inverse().matrix()) * vector_rep());
return HyperplaneHd<RT,LA>(dimension(),res.begin(),res.end()); }
/*{\Mtext \headerline{Non-Member Functions}}*/
static Comparison_result weak_cmp(
const HyperplaneHd<RT,LA>&, const HyperplaneHd<RT,LA>&);
@ Weak equality considers two hyperplanes equal if their coefficient
vectors are multiples of each other. We define the weak linear order
as the lexicographic order under weak equality. Let $i$ be minimal
such that either $h1_i$ or $h2_i$ is non-zero. We may assume that a
non-zero value is positive (since we consider weak equality). Thus if
exactly one of the vlaue is non-zero, we can decide the order right
there: The vector with the entry zero is smaller. If both entries are
non-zero, we compute scaling factors that make the $i$-th coefficients
equal and positive and proceed.
<<implementing HyperplaneHd>>=
template <class RT, class LA>
Comparison_result HyperplaneHd<RT,LA>::
weak_cmp(const HyperplaneHd<RT,LA>& h1,
const HyperplaneHd<RT,LA>& h2)
{
CGAL_assertion_msg((h1.dimension()==h2.dimension()),
"HyperplaneHd::weak_cmp: dimensions disagree.");
if(h1.identical(h2)) return EQUAL;
int i, d = h1.dimension();
for (i = 0; i <= d &&
h1.coefficient(i) == 0 &&
h2.coefficient(i) == 0; i++); // no body
if (h1.coefficient(i) == 0) return SMALLER;
if (h2.coefficient(i) == 0) return LARGER;
int s = CGAL_NTS sign(h1.coefficient(i)) *
CGAL_NTS sign(h2.coefficient(i));
RT s1 = (RT)s * h2.coefficient(i);
RT s2 = (RT)s * h1.coefficient(i);
// |s1 * h1.coefficient(i)| is
// $\Labs{ |h1.coefficient(i)*h2.coefficient(i)| }$
Comparison_result c;
while (++i <= d) {
c = CGAL_NTS compare(s1 * h1.coefficient(i),
s2 * h2.coefficient(i));
if (c != EQUAL) return c;
}
return EQUAL;
}
<<defining HyperplaneHd>>=
static Comparison_result strong_cmp(
const HyperplaneHd<RT,LA>&, const HyperplaneHd<RT,LA>&);
@ Strong equality considers two hyperplanes equal if their coefficient
vectors are positive multiples of each other. We define the strong
linear order as the lexicographic order under strong equality. Let $i$
be minimal such that either $h1_i$ or $h2_i$ is non-zero. If the
values have different signs we can decide the order right there: The
vector with the smaller entry is smaller. If the entries have the
same sign we compute positive scaling factors that make the $i$-th
coefficients equal and proceed.
<<implementing HyperplaneHd>>=
template <class RT, class LA>
Comparison_result HyperplaneHd<RT,LA>::
strong_cmp(const HyperplaneHd<RT,LA>& h1,
const HyperplaneHd<RT,LA>& h2)
{
CGAL_assertion_msg((h1.dimension()==h2.dimension()),
"HyperplaneHd::strong_cmp: dimensions disagree.");
if (h1.identical(h2)) return EQUAL;
int i;
int d = h1.dimension();
for (i = 0; i <=d &&
h1.coefficient(i)==0 &&
h2.coefficient(i)==0; i++) ; // no body
int c1 = CGAL_NTS sign(h1.coefficient(i));
int c2 = CGAL_NTS sign(h2.coefficient(i));
if (c1 != c2) return CGAL_NTS compare(c1,c2);
RT s1 = (RT)CGAL_NTS sign(h2.coefficient(i)) * h2.coefficient(i);
RT s2 = (RT)CGAL_NTS sign(h1.coefficient(i)) * h1.coefficient(i);
Comparison_result c;
while (++i <= d) {
c = CGAL_NTS compare(s1 * h1.coefficient(i),
s2 * h2.coefficient(i));
if (c != EQUAL) return c;
}
return EQUAL;
}
<<defining HyperplaneHd>>=
bool operator==(const HyperplaneHd<RT,LA>& h2) const
{ if (identical(h2)) return true;
if (dimension()!=h2.dimension()) return false;
return HyperplaneHd<RT,LA>::strong_cmp(*this,h2) == EQUAL;
}
bool operator!=(const HyperplaneHd<RT,LA>& h2) const
{ return !operator==(h2); }
friend std::istream& operator>> <>
(std::istream&, HyperplaneHd<RT,LA>&);
friend std::ostream& operator<< <>
(std::ostream&, const HyperplaneHd<RT,LA>&);
}; // end of class HyperplaneHd
template <class RT, class LA>
bool weak_equality(const HyperplaneHd<RT,LA>& h1,
const HyperplaneHd<RT,LA>& h2)
/*{\Mfunc test for weak equality. }*/
{ if (h1.identical(h2)) return true;
if (h1.dimension()!=h2.dimension()) return false;
return HyperplaneHd<RT,LA>::weak_cmp(h1,h2) == EQUAL;
}
<<prototyping>>=
template <class RT, class LA>
std::istream& operator>>(std::istream&, HyperplaneHd<RT,LA>&);
template <class RT, class LA>
std::ostream& operator<<(std::ostream&, const HyperplaneHd<RT,LA>&);
<<implementing HyperplaneHd>>=
template <class RT, class LA>
std::istream& operator>>(std::istream& I, HyperplaneHd<RT,LA>& h)
{ h.copy_on_write(); h.ptr->read(I); return I; }
template <class RT, class LA>
std::ostream& operator<<(std::ostream& O, const HyperplaneHd<RT,LA>& h)
{ h.ptr->print(O,"HyperplaneHd"); return O; }
template <class RT, class LA>
inline CGAL::io_Operator io_tag(const HyperplaneHd<RT,LA>&)
{ return CGAL::io_Operator(); }
<<defining HyperplaneHd>>=
/*{\Mimplementation
Hyperplanes are implemented by arrays of integers as an item type.
All operations like creation, initialization, tests, vector
arithmetic, input and output on a hyperplane $h$ take time
$O(|h.dimension()|)$. coordinate access and |dimension()| take
constant time. The space requirement is $O(|h.dimension()|)$. }*/
//----------------------- end of file ----------------------------------
@ \end{document}
Reference in New Issue View Git Blame Copy Permalink