\( \newcommand{\E}{\mathrm{E}} \) \( \newcommand{\A}{\mathrm{A}} \) \( \newcommand{\R}{\mathrm{R}} \) \( \newcommand{\N}{\mathrm{N}} \) \( \newcommand{\Q}{\mathrm{Q}} \) \( \newcommand{\Z}{\mathrm{Z}} \) \( \def\ccSum #1#2#3{ \sum_{#1}^{#2}{#3} } \def\ccProd #1#2#3{ \sum_{#1}^{#2}{#3} }\)
CGAL 4.13 - Bounding Volumes
Approximate_min_ellipsoid_d/ellipsoid_for_maple.cpp
// Usage: ./maple_example > maple.text
// Then enter in Maple 'read "maple.text";'.
#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 Traits::Point Point;
typedef std::vector<Point> Point_list;
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<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";
}
}