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 ddimensional Euclidean space ^{d}, d ≥ 2.
An ellipsoid in ^{d} is a Cartesian pointset of the form { x ∈ ^{d}  x^{T} E x + x^{T} e + η ≤ 0 }, where E is some positive definite matrix from the set ℝ^{d × d}, e is some real dvector, and η ∈ ℝ. A pointset P ⊆ ^{d} is called fulldimensional if its affine hull has dimension d. For a finite, fulldimensional 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 fulldimensional pointset P ⊂ ^{d} and a real number ε ≥ 0, we say that an ellipsoid E ⊂ ^{d} is an (1+ε)appoximation to (P) if P ⊂ E and (E) ≤ (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 fulldimensional 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 fulldimensional 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 fulldimensional.)
If P is not fulldimensional, 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) fulldimensional pointset; once S is determined, the algorithm can be invoked again to compute an approximation to (the lowerdimensional) (P) in S. Since is_full_dimensional() might (due to rounding errors, see above) return false even though P is fulldimensional, the lowerdimensional 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 minimumvolume 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.
The ellipsoid E computed by the algorithm satisfies the inclusions
 E ⊆ (P) ⊆ E 
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 ddimensional 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 3dimensional Cgal kernel, respectively.
Approximate_min_ellipsoid_d<Traits>::FT  
typedef Traits::FT FT (which is always a
typedef to double).
 
Approximate_min_ellipsoid_d<Traits>::ET  
typedef Traits::ET ET (which is an exact number type used for exact computation like for example in achieved_epsilon()).
 
Approximate_min_ellipsoid_d<Traits>::Point  
typedef Traits::Point Point
 
Approximate_min_ellipsoid_d<Traits>::Cartesian_const_iterator  
typedef Traits::Cartesian_const_iterator Cartesian_const_iterator
 
Approximate_min_ellipsoid_d<Traits>::Center_coordinate_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().
 
Approximate_min_ellipsoid_d<Traits>::Axes_lengths_iterator  
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().
 
Approximate_min_ellipsoid_d<Traits>::Axis_direction_iterator  
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.
template < class Iterator >  
Approximate_min_ellipsoid_d<Traits> ame ( double eps, Iterator first, Iterator last, Traits traits = Traits());  
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.

The following methods can be used to query the achieved approximation ratio 1+ε' and the computed ellipsoid E = { x ∈ ^{d}  x^{T} E x + x^{T} 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 E might not be exactly representable as double numbers.)
unsigned int  ame.number_of_points () const  returns the number of points of ame, i.e., P.  
double  ame.achieved_epsilon () const 
returns a number
ε' such that the computed approximation is (under exact
arithmetic) guaranteed to be an (1+ε')approximation to
(P).
 
double  ame.defining_matrix ( int i, int j) const  
gives access to the (i,j)th entry of the matrix E in the
representation { x ∈ ^{d}  x^{T} E x + x^{T} e + η ≤ 0
} of the computed approximation ellipsoid E. The number returned by
this routine is (1+ε')(d+1) E_{ij}, where ε' is
the number returned by achieved_epsilon().
 
double  ame.defining_vector ( int i) const  
gives access to the ith entry of the vector e in the
representation { x ∈ ^{d}  x^{T} E x + x^{T} e + η ≤ 0
} of the computed approximation ellipsoid E. The number returned by
this routine is (1+ε')(d+1) e_{i}, where ε' is
the number returned by achieved_epsilon().
 
double  ame.defining_scalar () const 
gives access to the scalar η from the
representation { x ∈ ^{d}  x^{T} E x + x^{T} e + η ≤ 0
} of the computed approximation ellipsoid E. The number returned by
this routine is (1+ε')(d+1) (η+1), where ε' is
the number returned by achieved_epsilon().
 
Traits  ame.traits () const  returns a const reference to the traits class object.  
int  ame.dimension () const  returns the dimension of the ambient space, i.e., the dimension of the points P. 
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.
Center_coordinate_iterator  ame.center_cartesian_begin () 
returns an iterator pointing to the first of the d Cartesian
coordinates of the computed ellipsoid's center. The returned point is a floatingpoint approximation to the ellipsoid's exact center; no guarantee is given w.r.t. the involved relative error.
 
Center_coordinate_iterator  ame.center_cartesian_end () 
returns the pasttheend iterator corresponding to
center_cartesian_begin().
 
Axes_lengths_iterator  ame.axes_lengths_begin () 
returns an iterator pointing to the first of the d descendantly
sorted lengths of the computed ellipsoid's axes. The d returned
numbers are floatingpoint approximations to the exact
axeslengths of the computed ellipsoid; no guarantee is given
w.r.t. the involved relative error. (See also method
axes_direction_cartesian_begin().)
 
Axes_lengths_iterator  ame.axes_lengths_end () 
returns the pasttheend iterator corresponding to
axes_lengths_begin().
 
Axes_direction_coordinate_iterator  
ame.axis_direction_cartesian_begin ( int i)  
returns an iterator
pointing to the first of the d Cartesian coordinates of the
computed ellipsoid's ith axis direction (i.e., unit vector in
direction of the ellipsoid's ith axis). The direction described
by this iterator is a floatingpoint approximation to the exact
axis direction of the computed ellipsoid; no guarantee is given
w.r.t. the involved relative error. An approximation to the
length of axis i is given by the ith entry of
axes_lengths_begin().
 
Axes_direction_coordinate_iterator  
ame.axis_direction_cartesian_end ( int i)  
returns the pasttheend
iterator corresponding to
axis_direction_cartesian_begin().

An object ame is valid iff
We implement Khachyian's algorithm for rounding polytopes [Kha96]. Internally, we use doublearithmetic and (initially a single) Choleskydecomposition. The algorithm's running time is O(nd^{2}(ε^{1}+ln d + ln ln (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.01approximation of (P) as an EPSfile, which you can view using gv, for instance. (In both examples you can change the variables n and d to experiment with the code.)
File: examples/Approximate_min_ellipsoid_d/ellipsoid.cpp
#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> #include <vector> #include <iostream> 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_cube_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.
File: examples/Approximate_min_ellipsoid_d/ellipsoid_for_maple.cpp
#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> #include <vector> #include <iostream> #include <iomanip> 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_cube_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<d1) cout << ","; } cout << "]"; if (i<d1) cout << ","; cout << "\n"; } cout << "]);\n" << "mp:= Vector(["; for (int i=0; i<d; ++i) { cout << mel.defining_vector(i)/alpha; if (i<d1) cout << ","; } cout << "]);\n" << "eta:= " << (mel.defining_scalar()/alpha1.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. 