\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$. <>= template 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>= template 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. <>= template 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. <>= template struct Center_of_sphereCd { typedef typename R::Point_d Point_d; typedef typename R::FT FT; typedef typename R::LA LA; template 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: <>= template 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. <>= template 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|. 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. <>= template struct Barycentric_coordinatesCd { typedef typename R::Point_d Point_d; typedef typename R::LA LA; typedef typename R::FT FT; template 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|. 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. <>= template struct OrientationCd { typedef typename R::Point_d Point_d; typedef typename R::LA LA; template 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|. 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). <>= template struct Side_of_oriented_sphereCd { typedef typename R::Point_d Point_d; typedef typename R::LA LA; typedef typename R::FT FT; template 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| and $|orientation(first,last)| \neq |ZERO|$. <>= template struct Side_of_bounded_sphereCd { typedef typename R::Point_d Point_d; template 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| 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. <>= template struct Contained_in_simplexCd { typedef typename R::Point_d Point_d; typedef typename R::LA LA; typedef typename R::FT FT; template 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|. 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$. <>= template struct Contained_in_affine_hullCd { typedef typename R::Point_d Point_d; typedef typename R::LA LA; template 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|. 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$. <>= template struct Affine_rankCd { typedef typename R::Point_d Point_d; typedef typename R::Vector_d Vector_d; typedef typename R::LA LA; template 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|. A set of points is affinely independent if their affine rank is equal to the number of points minus 1. <>= template struct Affinely_independentCd { typedef typename R::Point_d Point_d; typedef typename R::LA LA; template 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. <>= template 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|. 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. <>= template struct Contained_in_linear_hullCd { typedef typename R::LA LA; typedef typename R::FT FT; typedef typename R::Vector_d Vector_d; template 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|. 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. <>= template struct Linear_rankCd { typedef typename R::LA LA; template 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|. A set of vectors is linearly independent if their linear rank is equal to the number of vectors in the set. <>= template struct Linearly_independentCd { typedef typename R::LA LA; template 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|. To find a basis within |A| we use the corresponding function |independent_columns| of our linear algebra package. <>= template struct Linear_baseCd { typedef typename R::LA LA; typedef typename R::FT FT; typedef typename R::Vector_d Vector_d; template 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 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; } }; <>= // ====================================================================== // // 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 // 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 #include #include #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 namespace CGAL { <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> } //namespace CGAL #endif //CGAL_FUNCTION_OBJECTSCD_H @ \end{document}