21 template<
typename Obj, 
typename OTYPE>
 
   22 tlp::Circle<Obj, OTYPE> & tlp::Circle<Obj, OTYPE>::merge(
const tlp::Circle<Obj, OTYPE> &c) {
 
   23   Vector<Obj,2, OTYPE> p1(*
this);
 
   24   Vector<Obj,2, OTYPE> p2(c);
 
   25   Vector<Obj,2, OTYPE> c12=p2-p1;
 
   26   double norm=c12.norm();
 
   28   if (norm < 0.0000001) {
 
   35   double dx=c12[0]/norm;
 
   36   double dy=c12[1]/norm;
 
   41   Obj tmpDist=p1.dist(p2)/2.;
 
   43   if ( (tmpDist<radius) || (tmpDist<c.radius)) {
 
   59 template<
typename Obj, 
typename OTYPE>
 
   60 bool tlp::Circle<Obj, OTYPE>::isIncludeIn(
const tlp::Circle<Obj, OTYPE> &c)
 const {
 
   61   Vector<Obj,2, OTYPE> dir=c-*
this;
 
   62   return (dir.norm()+radius)<=c.radius;
 
   65 template<
typename Obj, 
typename OTYPE>
 
   66 tlp::Circle<Obj, OTYPE> tlp::enclosingCircle(
const tlp::Circle<Obj, OTYPE> &c1,
const tlp::Circle<Obj, OTYPE> &c2) {
 
   67   Vector<Obj,2, OTYPE> dir=c2-c1;
 
   71     return Circle<Obj, OTYPE>(c1, std::max(c1.radius,  c2.radius));
 
   74   Vector<Obj, 2, OTYPE> ext1=c1-dir*c1.radius;
 
   75   Vector<Obj, 2, OTYPE> ext2=c2+dir*c2.radius;
 
   76   return Circle<Obj, OTYPE>((ext1+ext2)/Obj(2),(ext2-ext1).norm()/Obj(2));
 
   79 template<
typename Obj, 
typename OTYPE>
 
   80 tlp::Circle<Obj, OTYPE> tlp::lazyEnclosingCircle(
const std::vector<tlp::Circle<Obj, OTYPE> > & circles) {
 
   82   tlp::Vector<Obj,4, OTYPE> boundingBox;
 
   84   typename std::vector< tlp::Circle<Obj, OTYPE> >::const_iterator it=circles.begin();
 
   85   boundingBox[0]=(*it)[0]-(*it).radius;
 
   86   boundingBox[1]=(*it)[1]-(*it).radius;
 
   87   boundingBox[2]=(*it)[0]+(*it).radius;
 
   88   boundingBox[3]=(*it)[1]+(*it).radius;
 
   91   for (; it!=circles.end(); ++it) {
 
   92     boundingBox[0] = std::min(boundingBox[0] , ((*it)[0]-(*it).radius));
 
   93     boundingBox[1] = std::min(boundingBox[1] , ((*it)[1]-(*it).radius));
 
   94     boundingBox[2] = std::max(boundingBox[2] , ((*it)[0]+(*it).radius));
 
   95     boundingBox[3] = std::max(boundingBox[3] , ((*it)[1]+(*it).radius));
 
   98   tlp::Vector<Obj,2, OTYPE> center;
 
   99   center[0]=(boundingBox[0]+boundingBox[2])/2.;
 
  100   center[1]=(boundingBox[1]+boundingBox[3])/2.;
 
  101   Obj radius = std::max( (boundingBox[2]-boundingBox[0])/2., (boundingBox[3]-boundingBox[1])/2. );
 
  102   tlp::Circle<Obj, OTYPE> result(center,radius);
 
  105   for (
typename std::vector< tlp::Circle<Obj, OTYPE> >::const_iterator it=circles.begin(); it!=circles.end(); ++it)
 
  111 template<
typename Obj, 
typename OTYPE>
 
  112 tlp::Circle<Obj, OTYPE> tlp::enclosingCircle(
const std::vector<tlp::Circle<Obj, OTYPE> > & circles) {
 
  113   class OptimumCircleHull {
 
  114     const std::vector<tlp::Circle<Obj, OTYPE> > *circles;
 
  115     std::vector<unsigned> enclosedCircles;
 
  118     tlp::Circle<Obj, OTYPE> result;
 
  120     static tlp::Circle<Obj, OTYPE> enclosingCircle(
const tlp::Circle<Obj, OTYPE> &c1,
const tlp::Circle<Obj, OTYPE> &c2,
const tlp::Circle<Obj, OTYPE> &c3) {
 
  130       Obj b = -b1*r2*a1*a1*b3+a3*a3*r2*r2*r2+b2*b2*r3*r3*r3+b1*b1*r3*r3*r3+b1*b1*r2*r2*r2+r1*r1*r1*b3*b3+r1*r1*r1*b2*b2+r2*r2*r2*b3*b3
 
  131               +a1*a1*r3*r3*r3+a2*a2*r1*r1*r1+a2*a2*r3*r3*r3+a1*a1*r2*r2*r2+a3*a3*r1*r1*r1-b2*b2*r3*a3*a3-b2*b2*r3*b3*b3-2*b2*r3*r3*r3*b1
 
  132               -b1*b1*r3*a3*a3-b1*b1*r3*b3*b3-b1*b1*r2*a2*a2+2*b1*b1*r2*b3*b3-b1*b1*r2*r3*r3-r1*b3*b3*b3*b2+r1*b3*b3*b3*b1
 
  133               +2*r1*b2*b2*b3*b3-r1*b2*b2*r3*r3+r2*b3*b3*b3*b2-r2*b3*b3*b3*b1-b2*b2*b2*r3*b1+2*b2*b2*r3*b1*b1+b2*b2*b2*r3*b3
 
  134               -b2*b2*r3*r1*r1+b1*b1*b1*r3*b3-b1*b1*b1*r3*b2-b1*b1*r3*r2*r2-b1*b1*r2*b2*b2-b1*b1*b1*r2*b3+b1*b1*b1*r2*b2
 
  135               -2*b1*r2*r2*r2*b3-r1*b3*b3*b1*b1-2*r1*r1*r1*b3*b2-r1*b3*b3*r2*r2-r1*b3*b3*a1*a1+r1*b2*b2*b2*b1
 
  136               -r1*b2*b2*b1*b1-r1*b2*b2*b2*b3-r1*b2*b2*a1*a1-r2*b3*b3*a2*a2-r2*b3*b3*b2*b2-r2*b3*b3*r1*r1
 
  137               -b2*r3*b1*b1*b3+b2*r3*a2*a2*b3+b2*r3*r1*r1*b3-b2*r3*r2*r2*b3+b2*r3*a1*a1*b3-2*b2*r3*a1*a2*b3+a1*a3*b2*b2*r3
 
  138               -2*a1*a3*b2*b1*r3-2*a1*a3*b2*b1*r2-2*a1*a3*b2*r1*b3+a1*a3*b2*b2*r1-2*a1*a3*b2*r2*b3-b2*r3*b1*a2*a2
 
  139               +2*b2*r3*b1*a3*a3+2*b2*r3*b1*b3*b3+b1*r2*b2*a3*a3-b1*r2*b2*b3*b3+b1*r2*b2*r3*r3+r1*b3*b1*a2*a2
 
  140               -r1*b3*b2*a3*a3+r1*b3*b2*r3*r3+r1*b3*b1*a3*a3-r1*b3*b1*r3*r3+r1*b2*b1*a2*a2+r1*b2*b1*a3*a3
 
  141               -r1*b2*b1*b3*b3+r1*b2*b1*r3*r3+2*r2*b3*b1*a2*a2+r2*b3*b2*a3*a3-r2*b3*b2*r3*r3-r2*b3*b1*a3*a3
 
  142               +r2*b3*b1*r3*r3+b2*r3*b1*r2*r2+4*b2*r3*a1*a2*b1+b1*r3*a2*a2*b3-b1*r3*b2*b2*b3
 
  143               -b1*r3*r1*r1*b3+b1*r3*r1*r1*b2+b1*r3*r2*r2*b3+b1*r3*a1*a1*b3-b1*r3*a1*a1*b2-2*b1*r3*a1*a2*b3-b1*b1*r3*a1*a2
 
  144               +b1*b1*r3*a1*a3+2*b1*r2*b2*b2*b3+b1*r2*r1*r1*b3-b1*r2*r1*r1*b2+b1*r2*a1*a1*b2
 
  145               -2*b1*r2*a1*a2*b3+b1*b1*r2*a1*a2-b1*b1*r2*a1*a3-r1*b3*b1*b2*b2+2*r1*b3*b1*b1*b2+2*r1*b3*a1*a1*b2+r1*b3*b1*r2*r2
 
  146               +r1*b3*b3*a1*a2-r1*b2*a2*a2*b3+r1*b2*r2*r2*b3-r1*b2*b1*r2*r2-2*r1*b2*a1*a2*b3-r2*b3*b1*b1*b2
 
  147               +r2*b3*r1*r1*b2+r2*b3*a1*a1*b2+r2*b3*b3*a1*a2+4*r2*b3*a1*a3*b1-a1*a1*a3*a3*r3+2*a1*a1*a3*a3*r2+a1*a3*a3*a3*r1
 
  148               -a1*a1*b3*b3*r3-a1*b3*b3*a3*r2+a1*r3*r3*a3*r2-a2*a1*a1*a3*r2+a2*b1*b1*a3*r2+a2*b3*b3*a3*r2
 
  149               -a1*a3*a3*a2*r1+2*a1*a3*a3*a2*r3+2*a1*b3*b3*a2*r3+a1*b3*b3*a3*r1-2*a1*r3*r3*r3*a2-a1*a1*r3*r3*r2
 
  150               -a2*a1*a1*a1*r3-a2*a2*a1*a1*r1+2*a2*a2*a1*a1*r3+a2*a1*a1*a1*r2-a2*a2*b1*b1*r1+2*a2*a2*a3*a3*r1-a2*a2*a3*a3*r3
 
  151               -a2*a3*a3*a3*r1-a2*a2*b3*b3*r3-a2*a2*r1*r1*r3-2*a2*r1*r1*r1*a3+a1*r3*r3*a2*r1-a1*r3*r3*a3*r1
 
  152               +2*a2*a1*a1*a3*r1+2*a2*b1*b1*a3*r1-a2*a3*a3*a1*r2-a2*b3*b3*a3*r1+a2*r1*r1*a1*r3-a2*r1*r1*a1*r2-a2*a2*r3*r3*r1
 
  153               -a1*a3*a3*a3*r2+a2*a3*a3*a3*r2-a1*a1*r3*r2*r2+a1*a1*a1*r3*a3+a2*a2*a2*r1*a1-a2*a2*a2*r1*a3-a2*a2*a2*r3*a1+a2*a2*a2*r3*a3
 
  154               -a1*a1*r2*a2*a2-a1*a1*r2*b2*b2-a1*a1*a1*r2*a3-2*a1*r2*r2*r2*a3-a3*a3*r1*a1*a1-a3*a3*r1*b1*b1-a3*a3*r1*r2*r2
 
  155               -a3*a3*r2*a2*a2-a3*a3*r2*b2*b2-a3*a3*r2*r1*r1+a2*r3*r3*a1*r2+a2*r3*r3*a3*r1-2*a3*r2*b1*a2*b3
 
  156               +a2*r1*r1*a3*r2-a2*r3*r3*a3*r2-a1*r3*a3*a2*a2-a1*r3*a3*r1*r1+a1*r3*a3*r2*r2+a2*r1*a1*b2*b2
 
  157               -a2*r1*a1*r2*r2-a2*r1*a3*b2*b2+a2*r1*a3*r2*r2-2*a2*r1*b1*a3*b2-a2*r3*a1*b2*b2+a2*r3*a1*r2*r2
 
  158               -a2*r3*a3*a1*a1+a2*r3*a3*b1*b1+a2*r3*a3*b2*b2+a2*r3*a3*r1*r1-a2*r3*a3*r2*r2-2*a2*r3*b1*a3*b2+2*a1*r2*a3*a2*a2
 
  159               +2*a1*r2*a3*b2*b2+a1*r2*a3*r1*r1-a3*r1*a1*a2*a2+a3*r1*a1*r2*r2-2*a3*r1*b1*a2*b3
 
  161       Obj tmp=a2*b3-a3*b2-a2*b1-a1*b3+a1*b2+a3*b1;
 
  162       Obj d = sqrt((a2*a2+b2*b2-r2*r2-2*a2*a3+a3*a3-r3*r3-2*b3*b2+b3*b3+2*r3*r2)*
 
  163                    (b3*b3+b1*b1-r1*r1-r3*r3-2*b1*b3+2*r3*r1+a3*a3-2*a1*a3+a1*a1)*
 
  164                    (a1*a1-2*a1*a2-r1*r1+b1*b1+b2*b2-r2*r2-2*b2*b1+2*r2*r1+a2*a2)*
 
  168       if (v<0) 
return tlp::Circle<Obj, OTYPE>(0,0,0);
 
  170       Obj aa = -2*a3*b2*b2*a1-2*a2*b3*b3*a1
 
  171                -2*a1*a1*b3*b2+a3*a3*b2*b2+a2*a2*b3*b3-r1*r1*b3*b3-r1*r1*b2*b2+a2*a2*b1*b1+a3*a3*b1*b1-a2*a2*r1*r1
 
  172                -a2*a2*r3*r3-a3*a3*r2*r2-2*a3*b2*a2*b3+2*a3*b2*a2*b1+2*a2*b3*a3*b1-2*b1*a2*a2*b3-a3*a3*r1*r1
 
  173                -2*a3*a3*b2*b1+2*b1*r3*r3*b2+2*r1*r1*b3*b2+2*b3*b1*r2*r2-2*a2*b1*b1*a3+2*a2*a3*r1*r1-b2*b2*r3*r3-b1*b1*r3*r3
 
  174                -b1*b1*r2*r2-r2*r2*b3*b3-a1*a1*r3*r3-a1*a1*r2*r2+a1*a1*b3*b3+a1*a1*b2*b2+2*b2*b2*r3*r1+2*b1*b1*r3*r2
 
  175                +2*r1*b3*b3*r2-2*b2*r3*b1*r2-2*b2*r3*r1*b3+2*b2*r3*r2*b3+2*b1*r3*r1*b3-2*b1*r3*r1*b2
 
  176                -2*b1*r3*r2*b3-2*b1*r2*r1*b3+2*b1*r2*r1*b2-2*r1*b2*r2*b3+2*a1*r3*r3*a2+2*a1*a1*r3*r2+2*a2*a2*r1*r3
 
  177                +2*a1*r2*r2*a3+2*a3*a3*r1*r2-2*a1*r3*a2*r1+2*a1*r3*a3*r1+2*a2*r1*a1*r2-2*a2*r3*a1*r2
 
  178                -2*a2*r3*a3*r1-2*a1*r2*a3*r1-2*a1*r3*a3*r2-2*a2*r1*a3*r2+2*a2*r3*a3*r2+2*a2*b3*a1*b2
 
  179                +2*a3*b2*a1*b3+2*a2*b1*a1*b3-2*a2*b1*a1*b2-2*a1*b3*a3*b1+2*a1*b2*a3*b1;
 
  181       Obj y = -0.5*(-2*a1*R*r2-2*a3*R*r1+2*a2*R*r1+2*a3*R*r2-2*a2*R*r3+a1*a3*a3+a1*b3*b3-a1*r3*r3
 
  182                     +a2*a1*a1+a2*b1*b1+a1*r2*r2+2*a1*R*r3-a3*b1*b1+a3*a2*a2+a3*b2*b2+a3*r1*r1-a3*r2*r2-a2*a3*a3-a2*b3*b3
 
  183                     -a2*r1*r1+a2*r3*r3-a1*a2*a2-a1*b2*b2-a3*a1*a1)/(a2*b3+a3*b1-a2*b1-a1*b3+a1*b2-a3*b2);
 
  184       Obj x = 0.5*(-a1*a1*b3+a1*a1*b2+2*R*r2*b3+b1*a3*a3+b1*b3*b3+2*b1*R*r3-2*R*r1*b3+2*R*r1*b2+b1*r2*r2
 
  185                    -b2*a3*a3-b2*b3*b3+b2*r3*r3-2*b2*R*r3-b1*r3*r3-r2*r2*b3+a2*a2*b3-r1*r1*b2-b1*b2*b2+b1*b1*b2-b1*b1*b3
 
  186                    -b1*a2*a2+b2*b2*b3+r1*r1*b3-2*b1*R*r2)/(a2*b3+a3*b1-a2*b1-a1*b3+a1*b2-a3*b2);
 
  187       return  tlp::Circle<Obj, OTYPE>(x,y,R);
 
  192         result=tlp::enclosingCircle((*circles)[b1],(*circles)[b2]);
 
  195         unsigned selectedCircle=popBack();
 
  198         if (!(*circles)[selectedCircle].isIncludeIn(result)) {
 
  199           result=enclosingCircle((*circles)[b1],(*circles)[b2],(*circles)[selectedCircle]);
 
  200           pushFront(selectedCircle);
 
  203           pushBack(selectedCircle);
 
  210         result=(*circles)[b1];
 
  213         unsigned selectedCircle=popBack();
 
  216         if (!(*circles)[selectedCircle].isIncludeIn(result)) {
 
  219           pushFront(selectedCircle);
 
  222           pushBack(selectedCircle);
 
  228         result=Circle<Obj, OTYPE>(0,0,0);
 
  231         unsigned selectedCircle=popBack();
 
  234         if (!(*circles)[selectedCircle].isIncludeIn(result)) {
 
  237           pushFront(selectedCircle);
 
  240           pushBack(selectedCircle);
 
  244     bool isEmpty()
 const {
 
  245       return first==(last+1)%enclosedCircles.size();
 
  247     void pushFront(
unsigned c) {
 
  248       first=(first+enclosedCircles.size()-1)%enclosedCircles.size();
 
  249       enclosedCircles[first]=c;
 
  252       unsigned result=enclosedCircles[last];
 
  253       last=(last+enclosedCircles.size()-1)%enclosedCircles.size();
 
  256     void pushBack(
unsigned c) {
 
  257       last=(last+1)%enclosedCircles.size();
 
  258       enclosedCircles[last]=c;
 
  261     OptimumCircleHull():circles(NULL),first(0),last(0),b1(0),b2(0) {}
 
  263     tlp::Circle<Obj, OTYPE> operator()(
const std::vector<tlp::Circle<Obj, OTYPE> > &circlesSet) {
 
  265       enclosedCircles.resize(circlesSet.size()+1);
 
  267       last=circlesSet.size()-1;
 
  269       for (
unsigned i=0; i<circlesSet.size(); ++i)
 
  270         enclosedCircles[i]=i;
 
  272       for (
unsigned i=circlesSet.size(); i>0;) {
 
  273         unsigned idx=(unsigned)(1.0*rand()*i/RAND_MAX);
 
  275         std::swap(enclosedCircles[idx],enclosedCircles[i]);
 
  282   return OptimumCircleHull()(circles);
 
  285 template <
typename Obj, 
typename OTYPE>
 
  286 std::ostream& tlp::operator<<(std::ostream &os,const tlp::Circle<Obj, OTYPE> &a) {
 
  287   os << 
"((" << a[0] << 
"," <<  a[1] << 
")," << a.radius << 
")";