|
An object of class Approximate_min_ellipsoid_d<Traits> is an approximation to the
ellipsoid of smallest volume enclosing a finite multiset of points
in d-dimensional Euclidean space 
d, d 
2.
An ellipsoid in d is a Cartesian pointset of the form {
x 
d 
xT E x + xT e + 
0 }, where E is some
positive definite matrix from the set 
d × d, e is some
real d-vector, and . A pointset P 
d is
called full-dimensional if its affine hull has dimension d.
For a finite, full-dimensional pointset P we denote by (P) the
smallest ellipsoid that contains all points of P; this ellipsoid
exists and is unique.
For a given finite and full-dimensional pointset P 
d and a
real number 
0, we say that an ellipsoid 
d is an (1+)-appoximation to (P) if
P and () (1+)
((P)). In other words, an (1+)-approximation to
(P) is an enclosing ellipsoid whose volume is by at most a
factor of 1+ larger than the volume of the smallest
enclosing ellipsoid of P.
Given this notation, an object of class Approximate_min_ellipsoid_d<Traits> represents an
(1+)-approximation to (P) for a given finite and
full-dimensional multiset of points P d and a real constant
>0.1 When an
Approximate_min_ellipsoid_d<Traits> object is constructed, an
iterator over the points P and the number have to be
specified; the number defines the desired
approximation ratio 1+. The underlying algorithm will then
try to compute an (1+)-approximation to (P), and one of
the following two cases takes place.
Important note: due to rounding errors, the algorithm cannot
in all cases decide correctly whether P is full-dimensional or
not. If is_full_dimensional() returns false, the points
lie in such a ``thin'' subspace of d that the algorithm is
incapable of computing an approximation to (P). More
precisely, if is_full_dimensional() returns false, there
exist two parallel hyperplanes in d with the points P in
between so that the distance 
between the hyperplanes is
very small, possible zero. (If =0 then P is not
full-dimensional.)
If P is not full-dimensional, linear algebra techniques should be
used to determine an affine subspace S of d that contains the
points P as a (w.r.t. S) full-dimensional pointset; once S is
determined, the algorithm can be invoked again to compute an
approximation to (the lower-dimensional) (P) in S. Since
is_full_dimensional() might (due to rounding errors, see
above) return false even though P is full-dimensional, the
lower-dimensional subspace S containing P need not exist.
Therefore, it might be more advisable to fit a hyperplane H
through the pointset P, project P onto this affine subspace H,
and compute an approximation to the minimum-volume enclosing
ellipsoid of the projected points within H; the fitting can be
done for instance using the linear_least_squares_fitting()
function from the CGAL package Principal_component_analysis.
to (P), but
depending on the input problem (i.e., on the pair (P,)),
it may not have achieved the desired approximation ratio but merely
some worse approximation ratio 1+'>1+. The
achieved approximation ratio 1+' can be queried using
achieved_epsilon(), which returns '. The ellipsoid
itself can be queried via the methods
defining_matrix(), defining_vector(), and
defining_scalar().
The ellipsoid computed by the algorithm satisfies the inclusions
where f denotes the ellipsoid scaled by the
factor f + with respect to its center, and where (A) denotes the convex hull of a pointset
A d.
The underlying algorithm can cope with all kinds of inputs (multisets
P, [0, 
)) and terminates in all cases. There
is, however, no guarantee that any desired approximation ratio
is actually achieved; the performance of the algorithm in this respect
highly depends on the input pointset. Values of at least 0.01 for
are usually handled without problems.
Internally, the algorithm represents the input points' Cartesian
coordinates as double's. For this conversion to work, the input
point coordinates must be convertible to double. Also, in order
to compute the achieved epsilon ' mentioned above, the algorithm
requires a number type ET that provides exact arithmetic.
(Both these aspects are discussed in the documentation of the concept
ApproximateMinEllipsoid_d_Traits_d.)
#include <CGAL/Approximate_min_ellipsoid_d.h>
The template parameter Traits is a model for ApproximateMinEllipsoid_d_Traits_d.
We provide the model CGAL::Approximate_min_ellipsoid_d_traits_d<K> using the d-dimensional CGAL kernel; the models CGAL::Approximate_min_ellipsoid_d_traits_2<K> and CGAL::Approximate_min_ellipsoid_d_traits_3<K> are for use with the 2- and 3-dimensional CGAL kernel, respectively.
|
| |
|
typedef Traits::FT FT (which is always a
typedef to double).
| |
|
| |
|
typedef Traits::ET ET (which is an exact number type used for exact computation like for example in achieved_epsilon()).
| |
|
| |
|
typedef Traits::Point Point
| |
|
| |
|
typedef Traits::Cartesian_const_iterator Cartesian_const_iterator
| |
|
| |
|
A model of STL concept
RandomAccessIterator with value type double that is used
to iterate over the Cartesian center coordinates of the computed
ellipsoid, see center_cartesian_begin().
| |
|
| |
|
A model of STL concept
RandomAccessIterator with value type double that is used
to iterate over the lengths of the semiaxes of the computed ellipsoid,
see axes_lengths_begin().
| |
|
| |
|
A model of STL concept
RandomAccessIterator with value type double that is used
to iterate over the Cartesian coordinates of the direction of a fixed
axis of the computed ellipsoid, see
axis_direction_cartesian_begin().
| |
An object of type Approximate_min_ellipsoid_d<Traits> can be created from an arbitrary point set P and some nonnegative double value eps.
|
| |||
| |||
|
initializes ame to an (1+)-approximation of
(P) with P being the set of points in the range
[first,last). The number in this will
be at most eps, if possible. However, due to the
limited precision in the algorithm's underlying arithmetic, it
can happen that the computed approximation ellipsoid has a
worse approximation ratio (and can thus be larger
than eps in general). In any case, the number
(and with this, the achived approximation
1+) can be queried by calling the routine
achieved_epsilon() discussed below. Requirement: Iterator must be a model for concept InputIterator with value type Point. Precondition: The dimension d of the input points must be at least 2, and >0.
| |||
The following methods can be used to query the achieved approximation
ratio 1+' and the computed ellipsoid = { x d
xT E x + xT e + 0 }. The methods
defining_matrix(), defining_vector(), and
defining_scalar() do not return E, e, and directly
but yield multiples of these quantities that are exactly representable
using the double type. (This is necessary because the
parameters E, e, and of the computed approximation
ellipsoid might not be exactly representable as
double numbers.)
In order to access the center and semiaxes of the computed approximation ellipsoid, the functions center_cartesian_begin(), axes_lengths_begin(), and axis_direction_cartesian_begin() can be used. In constrast to the above access functions achieved_epsilon(), defining_matrix(), defining_vector(), and defining_scalar(), which return the described quantities exactly, the routines below return numerical approximations to the real center and real semiaxes of the computed ellipsoid; the comprised relative error may be larger than zero, and there are no guarantees for the returned quantities.
An object ame is valid iff
')-approximation to the smallest
ellipsoid (P) of P,
|
|
| |
| returns true iff ame is valid according to the above definition. If verbose is true, some messages concerning the performed checks are written to the standard error stream. | ||
We implement Khachyian's algorithm for rounding
polytopes [Kha96]. Internally, we use
double-arithmetic and (initially a single)
Cholesky-decomposition. The algorithm's running time is 
(nd2(-1+lnd + lnln(n))), where n=|P| and
1+ is the desired approximation ratio.
To illustrate the usage of Approximate_min_ellipsoid_d<Traits> we give two examples in 2D. The
first program generates a random set P 2 and outputs the
points and a 1.01-approximation of (P) as an EPS-file, which
you can view using gv, for instance. (In both examples you can
change the variables n and d to experiment with the code.)
#include <vector>
#include <iostream>
#include <CGAL/basic.h>
#include <CGAL/Cartesian_d.h>
#include <CGAL/MP_Float.h>
#include <CGAL/point_generators_d.h>
#include <CGAL/Approximate_min_ellipsoid_d.h>
#include <CGAL/Approximate_min_ellipsoid_d_traits_d.h>
typedef CGAL::Cartesian_d<double> Kernel;
typedef CGAL::MP_Float ET;
typedef CGAL::Approximate_min_ellipsoid_d_traits_d<Kernel, ET> Traits;
typedef Traits::Point Point;
typedef std::vector<Point> Point_list;
typedef CGAL::Approximate_min_ellipsoid_d<Traits> AME;
int main()
{
const int n = 1000; // number of points
const int d = 2; // dimension
const double eps = 0.01; // approximation ratio is (1+eps)
// create a set of random points:
Point_list P;
CGAL::Random_points_in_iso_box_d<Point> rpg(d,100.0);
for (int i = 0; i < n; ++i) {
P.push_back(*rpg);
++rpg;
}
// compute approximation:
Traits traits;
AME ame(eps, P.begin(), P.end(), traits);
// write EPS file:
if (ame.is_full_dimensional() && d == 2)
ame.write_eps("example.eps");
// output center coordinates:
std::cout << "Cartesian center coordinates: ";
for (AME::Center_coordinate_iterator c_it = ame.center_cartesian_begin();
c_it != ame.center_cartesian_end();
++c_it)
std::cout << *c_it << ' ';
std::cout << ".\n";
if (d == 2 || d == 3) {
// output axes:
AME::Axes_lengths_iterator axes = ame.axes_lengths_begin();
for (int i = 0; i < d; ++i) {
std::cout << "Semiaxis " << i << " has length " << *axes++ << "\n"
<< "and Cartesian coordinates ";
for (AME::Axes_direction_coordinate_iterator
d_it = ame.axis_direction_cartesian_begin(i);
d_it != ame.axis_direction_cartesian_end(i); ++d_it)
std::cout << *d_it << ' ';
std::cout << ".\n";
}
}
}
The second program outputs the approximation in a format suitable for display in Maplesoft's Maple.
// Usage: ./maple_example > maple.text
// Then enter in Maple 'read "maple.text";'.
#include <vector>
#include <iostream>
#include <iomanip>
#include <CGAL/basic.h>
#include <CGAL/Cartesian_d.h>
#include <CGAL/MP_Float.h>
#include <CGAL/point_generators_d.h>
#include <CGAL/Approximate_min_ellipsoid_d.h>
#include <CGAL/Approximate_min_ellipsoid_d_traits_d.h>
typedef CGAL::Cartesian_d<double> Kernel;
typedef CGAL::MP_Float ET;
typedef CGAL::Approximate_min_ellipsoid_d_traits_d<Kernel, ET> Traits;
typedef Traits::Point Point;
typedef std::vector<Point> Point_list;
typedef CGAL::Approximate_min_ellipsoid_d<Traits> AME;
int main()
{
const int n = 100; // number of points
const int d = 2; // dimension
const double eps = 0.01; // approximation ratio is (1+eps)
// create a set of random points:
Point_list P;
CGAL::Random_points_in_iso_box_d<Point> rpg(d,1.0);
for (int i = 0; i < n; ++i) {
P.push_back(*rpg);
++rpg;
}
// compute approximation:
Traits traits;
AME mel(eps, P.begin(), P.end(), traits);
// output for Maple:
if (mel.is_full_dimensional() && d == 2) {
const double alpha = (1+mel.achieved_epsilon())*(d+1);
// output points:
using std::cout;
cout << "restart;\n"
<< "with(LinearAlgebra):\n"
<< "with(plottools):\n"
<< "n:= " << n << ":\n"
<< "P:= Matrix(" << d << "," << n << "):\n";
for (int i=0; i<n; ++i)
for (int j=0; j<d; ++j)
cout << "P[" << j+1 << "," << i+1 << "] := "
<< std::setiosflags(std::ios::scientific)
<< std::setprecision(20) << P[i][j] << ":\n";
cout << "\n";
// output defining equation:
cout << "Mp:= Matrix([\n";
for (int i=0; i<d; ++i) {
cout << " [";
for (int j=0; j<d; ++j) {
cout << mel.defining_matrix(i,j)/alpha;
if (j<d-1)
cout << ",";
}
cout << "]";
if (i<d-1)
cout << ",";
cout << "\n";
}
cout << "]);\n" << "mp:= Vector([";
for (int i=0; i<d; ++i) {
cout << mel.defining_vector(i)/alpha;
if (i<d-1)
cout << ",";
}
cout << "]);\n"
<< "eta:= " << (mel.defining_scalar()/alpha-1.0) << ";\n"
<< "v:= Vector([x,y]):\n"
<< "e:= Transpose(v).Mp.v+Transpose(v).mp+eta;\n"
<< "plots[display]({seq(point([P[1,i],P[2,i]]),i=1..n),\n"
<< " plots[implicitplot](e,x=-5..5,y=-5..5,numpoints=10000)},\n"
<< " scaling=CONSTRAINED);\n";
}
}
| 1 | A multiset is a set where elements may have multiplicity greater than 1. |