MUON/AliMUONContourMaker.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 // $Id$
  17
  18 ///
  19 /// Maker/merger of contours. Can create contour from a set polygons or
  20 /// merger a set of contours into a single one.
  21 ///
  22 /// This is based on (one of the) algorithm found in
  23 /// Diane L. Souvaine and Iliana Bjorling-Sachs,
  24 /// Proceedings of the IEEE, Vol. 80, No. 9, September 1992, p. 1449
  25 ///
  26 /// Note that besides the AliMUON prefix, nothing is really MUON specific
  27 /// in this class...
  28 ///
  29 /// \author Laurent Aphecetche, Subatech
  30 ///
  31
  32 #include "AliMUONContourMaker.h"
  33
  34 #include "AliCodeTimer.h"
  35 #include "AliLog.h"
  36 #include "AliMUONContour.h"
  37 #include "AliMUONPointWithRef.h"
  38 #include "AliMUONPolygon.h"
  39 #include "AliMUONSegment.h"
  40 #include "AliMUONSegmentTree.h"
  41 #include "Riostream.h"
  42 #include "TArrayD.h"
  43 #include "TMath.h"
  44 #include <cassert>
  45
  46 /// \cond CLASSIMP
  47 ClassImp(AliMUONContourMaker)
  48 /// \endcond
  49
  50 namespace
  51 {
  52   void PrintSegments(const TObjArray& array)
  53   {
  54     TIter next(&array);
  55     AliMUONSegment* s;
  56     Int_t i(0);
  57     while ( ( s = static_cast<AliMUONSegment*>(next()) ) )
  58     {
  59       cout << Form("i=%d %s",i,s->AsString()) << endl;
  60       ++i;
  61     }
  62   }
  63 }
  64
  65 //_____________________________________________________________________________
  66 AliMUONContourMaker::AliMUONContourMaker()
  67 {
  68 }
  69
  70 //_____________________________________________________________________________
  71 AliMUONContourMaker::~AliMUONContourMaker()
  72 {
  73 }
  74
  75 //_____________________________________________________________________________
  76 AliMUONContour*
  77 AliMUONContourMaker::CreateContour(const TObjArray& polygons, const char* name) const
  78 {
  79   /// Create the contour of the polygon array
  80   /// and get back the intermediate verticals and horizontal segments
  81   /// both arrays are arrays of AliMUONSegment objects.
  82
  83   AliCodeTimerAuto("");
  84
  85   if ( polygons.IsEmpty() ) return 0x0; // protection against user error...
  86
  87   // Sanity check : insure that all polygons are oriented counter-clockwise
  88   TIter next(&polygons);
  89   AliMUONPolygon* pol;
  90   while ( ( pol = static_cast<AliMUONPolygon*>(next()) ) )
  91   {
  92     if ( !pol->IsCounterClockwiseOriented() )
  93     {
  94       AliError(Form("Got a clockwise oriented polygon in CreateContour(%s). That's not OK !",name));
  95       StdoutToAliError(polygons.Print());
  96       return 0x0;
  97     }
  98   }
  99
 100   AliMUONContour* contour(0x0);
 101
 102   if ( polygons.GetLast() == 0 )
 103   {
 104     AliCodeTimerAuto("Trivial case");
 105     contour = new AliMUONContour(name);
 106     pol = static_cast<const AliMUONPolygon*>(polygons.First());
 107     contour->Add(*pol);
 108     contour->AssertOrientation();
 109     return contour;
 110   }
 111
 112   TObjArray polygonVerticalEdges; // arrray of AliMUONSegment objects
 113   polygonVerticalEdges.SetOwner(kTRUE);
 114   // get vertical edges of input polygons
 115   GetVerticalEdges(polygons,polygonVerticalEdges);
 116
 117   // sort them in ascending x order
 118   // if same x, insure that left edges are before right edges
 119   // within same x, order by increasing bottommost y (see AliMUONSegment::Compare method)
 120   polygonVerticalEdges.Sort();
 121
 122   if ( polygonVerticalEdges.GetLast()+1 < 2 )
 123   {
 124     AliError(Form("Got too few edges here for createContour %s",name));
 125     polygons.Print();
 126     TObject* o(0x0);
 127     o->Print();
 128   }
 129
 130   // Find the vertical edges of the merged contour. This is the meat of the algorithm...
 131   TObjArray contourVerticalEdges;
 132   contourVerticalEdges.SetOwner(kTRUE);
 133   Sweep(polygonVerticalEdges,contourVerticalEdges);
 134
 135   TObjArray horizontals;
 136   horizontals.SetOwner(kTRUE);
 137   VerticalToHorizontal(contourVerticalEdges,horizontals);
 138
 139   contour = FinalizeContour(contourVerticalEdges,horizontals);
 140
 141   if ( contour && name ) contour->SetName(name);
 142
 143   return contour;
 144 }
 145
 146 //_____________________________________________________________________________
 147 AliMUONContour*
 148 AliMUONContourMaker::FinalizeContour(const TObjArray& verticals,
 149                                      const TObjArray& horizontals) const
 150 {
 151   /// For a list of vertical and horizontal edges, we build the final
 152   /// contour object.
 153
 154   AliCodeTimerAuto("");
 155
 156   TObjArray all; // array of AliMUONSegment
 157   TObjArray inorder; // array of AliMUONSegment
 158
 159   all.SetOwner(kFALSE);
 160   inorder.SetOwner(kFALSE);
 161
 162   for ( Int_t i = 0; i <= verticals.GetLast(); ++i )
 163   {
 164     all.Add(verticals.UncheckedAt(i));
 165     all.Add(horizontals.UncheckedAt(i));
 166   }
 167
 168   Int_t i(0);
 169
 170   AliMUONContour* contour = new AliMUONContour;
 171
 172   int total(0);
 173
 174   while ( !all.IsEmpty() )
 175   {
 176     total++;
 177
 178     if ( total > 1000 )
 179     {
 180       AliError("Total 1000 reached !!!!");
 181       return 0x0;
 182     }
 183
 184     AliMUONSegment* si = static_cast<AliMUONSegment*>(all.UncheckedAt(i));
 185     inorder.Add(si);
 186     const AliMUONSegment* all0 = static_cast<const AliMUONSegment*>(all.First());
 187     if ( i != 0 && AliMUONSegment::AreEqual(si->EndX(),all0->StartX()) && AliMUONSegment::AreEqual(si->EndY(),all0->StartY()) )
 188     {
 189       Int_t n(-1);
 190
 191       AliMUONPolygon polygon(inorder.GetLast()+2);
 192
 193       // we got a cycle. Add it to the contour
 194       for ( Int_t j = 0; j <= inorder.GetLast(); ++j )
 195       {
 196         AliMUONSegment* s = static_cast<AliMUONSegment*>(inorder.UncheckedAt(j));
 197         polygon.SetVertex(++n,s->StartX(),s->StartY());
 198         all.Remove(s);
 199       }
 200
 201       all.Compress();
 202
 203       polygon.Close();
 204
 205       contour->Add(polygon);
 206
 207       if ( ! all.IsEmpty() )
 208       {
 209         i = 0;
 210         inorder.Clear();
 211       }
 212       continue;
 213     }
 214
 215     for ( Int_t j = 0; j <= all.GetLast(); ++j)
 216     {
 217       if ( j != i )
 218       {
 219         const AliMUONSegment* sj = static_cast<const AliMUONSegment*>(all.UncheckedAt(j));
 220         if ( AliMUONSegment::AreEqual(si->EndX(),sj->StartX()) && AliMUONSegment::AreEqual(si->EndY(),sj->StartY()))
 221         {
 222           i = j;
 223           break;
 224         }
 225       }
 226     }
 227   }
 228
 229   contour->AssertOrientation(kTRUE);
 230   return contour;
 231 }
 232
 233
 234 //_____________________________________________________________________________
 235 void
 236 AliMUONContourMaker::GetVerticalEdges(const TObjArray& polygons, TObjArray& polygonVerticalEdges) const
 237 {
 238   /// From an array of polygons, extract the list of vertical edges.
 239   /// Output array polygonVerticalEdges should be empty before calling.
 240
 241   AliCodeTimerAuto("");
 242
 243   for ( Int_t i = 0; i <= polygons.GetLast(); ++i )
 244   {
 245     const AliMUONPolygon* g = static_cast<const AliMUONPolygon*>(polygons.UncheckedAt(i));
 246     for ( Int_t j = 0; j < g->NumberOfVertices()-1; ++j )
 247     {
 248       if ( AliMUONSegment::AreEqual(g->X(j),g->X(j+1)) ) // segment is vertical
 249       {
 250         polygonVerticalEdges.Add(new AliMUONSegment(g->X(j),g->Y(j),g->X(j+1),g->Y(j+1)));
 251       }
 252     }
 253   }
 254 }
 255
 256
 257 //_____________________________________________________________________________
 258 void
 259 AliMUONContourMaker::GetYPositions(const TObjArray& polygonVerticalEdges,
 260                                    TArrayD& yPositions) const
 261 {
 262   /// Fill the array yPositions with the different y positions found in
 263   /// polygonVerticalEdges
 264
 265   AliCodeTimerAuto("");
 266
 267   Double_t* y = new Double_t[polygonVerticalEdges.GetSize()*2];
 268   Int_t n(0);
 269
 270   for ( Int_t i = 0; i < polygonVerticalEdges.GetLast(); ++i )
 271   {
 272     AliMUONSegment* s = static_cast<AliMUONSegment*>(polygonVerticalEdges.UncheckedAt(i));
 273     y[n] = s->StartY();
 274     y[n+1] = s->EndY();
 275     n += 2;
 276   }
 277   Int_t* ix = new Int_t[n+1];
 278
 279   TMath::Sort(n,y,ix,kFALSE);
 280
 281   yPositions.Set(n+1);
 282
 283   Int_t u(0);
 284   Double_t x(FLT_MAX);
 285
 286   for ( Int_t i = 0; i < n; ++i )
 287   {
 288     if ( y[ix[i]] != x )
 289     {
 290       yPositions[u] = y[ix[i]];
 291       x = y[ix[i]];
 292       ++u;
 293     }
 294   }
 295
 296   yPositions.Set(u);
 297
 298   delete[] ix;
 299   delete[] y;
 300
 301 }
 302
 303 //_____________________________________________________________________________
 304 AliMUONContour*
 305 AliMUONContourMaker::MergeContour(const TObjArray& contours, const char* name) const
 306 {
 307   /// Merge all the polygons of all contours into a single contour
 308
 309   AliCodeTimerAuto("");
 310
 311   TObjArray polygons;
 312   polygons.SetOwner(kTRUE);
 313
 314   TIter next(&contours);
 315   AliMUONContour* contour;
 316   while ( ( contour = static_cast<AliMUONContour*>(next()) ) )
 317   {
 318     const TObjArray* contourPolygons = contour->Polygons();
 319     TIter nextPol(contourPolygons);
 320     AliMUONPolygon* pol;
 321     while ( ( pol = static_cast<AliMUONPolygon*>(nextPol()) ) )
 322     {
 323       polygons.Add(new AliMUONPolygon(*pol));
 324     }
 325   }
 326
 327   if ( polygons.IsEmpty() ) return 0x0;
 328
 329   contour = CreateContour(polygons,name);
 330
 331   return contour;
 332 }
 333
 334 //_____________________________________________________________________________
 335 void
 336 AliMUONContourMaker::SortPoints(const TObjArray& polygonVerticalEdges,
 337                                 TObjArray& sortedPoints) const
 338 {
 339   /// Sort the point of the vertical edges in ascending order, first on ordinate,
 340   /// then on abcissa, and put them in output vector sortedPoints.
 341   /// Output array sortedPoints should be empty before calling this method.
 342
 343   AliCodeTimerAuto("");
 344
 345   for ( Int_t i = 0; i <= polygonVerticalEdges.GetLast(); ++i )
 346   {
 347     const AliMUONSegment* e = static_cast<const AliMUONSegment*>(polygonVerticalEdges.UncheckedAt(i));
 348     sortedPoints.Add(new AliMUONPointWithRef(e->StartX(),e->StartY(),i));
 349     sortedPoints.Add(new AliMUONPointWithRef(e->EndX(),e->EndY(),i));
 350     // note that we keep track of the original edge, which is used
 351     // later on to deduce orientation of horizontal edges.
 352   }
 353
 354   sortedPoints.Sort();
 355 }
 356
 357 //_____________________________________________________________________________
 358 void
 359 AliMUONContourMaker::Sweep(const TObjArray& polygonVerticalEdges,
 360                            TObjArray& contourVerticalEdges) const
 361 {
 362   /// This is the meat of the algorithm of the contour merging...
 363
 364   AliCodeTimerAuto("");
 365
 366   TArrayD yPositions;
 367   GetYPositions(polygonVerticalEdges,yPositions);
 368
 369   AliMUONSegmentTree segmentTree(yPositions);
 370
 371   for ( Int_t i = 0; i <= polygonVerticalEdges.GetLast(); ++i )
 372   {
 373     const AliMUONSegment* edge = static_cast<const AliMUONSegment*>(polygonVerticalEdges.UncheckedAt(i));
 374
 375     assert(edge!=0x0);
 376
 377     if ( edge->IsLeftEdge() )
 378     {
 379       segmentTree.Contribution(edge->Bottom(),edge->Top());
 380       segmentTree.InsertInterval(edge->Bottom(),edge->Top());
 381     }
 382     else
 383     {
 384       segmentTree.DeleteInterval(edge->Bottom(),edge->Top());
 385       segmentTree.Contribution(edge->Bottom(),edge->Top());
 386     }
 387
 388     AliMUONSegment e1(*edge);
 389
 390     if ( i < polygonVerticalEdges.GetLast() )
 391     {
 392       const AliMUONSegment* next = static_cast<const AliMUONSegment*>(polygonVerticalEdges.UncheckedAt(i+1));
 393       e1 = *next;
 394     }
 395
 396     if ( ( edge->IsLeftEdge() != e1.IsLeftEdge() ) ||
 397         ( !AliMUONSegment::AreEqual(edge->StartX(),e1.StartX() ) ) ||
 398         ( i == polygonVerticalEdges.GetLast() ) )
 399     {
 400       const TObjArray& stack = segmentTree.Stack();
 401
 402       double x = edge->StartX();
 403
 404       for ( Int_t j = 0; j <= stack.GetLast(); ++j )
 405       {
 406         AliMUONSegment* sj = static_cast<AliMUONSegment*>(stack.UncheckedAt(j));
 407         AliMUONSegment* s = new AliMUONSegment(x,sj->StartY(),x,sj->EndY());
 408
 409         if  (s->IsAPoint())
 410         {
 411           delete s;
 412           continue;
 413         }
 414
 415         if ( edge->IsLeftEdge() != s->IsLeftEdge() )
 416         {
 417           s->Set(x,sj->EndY(),x,sj->StartY());
 418         }
 419         contourVerticalEdges.Add(s);
 420       }
 421       segmentTree.ResetStack();
 422     }
 423   }
 424 }
 425
 426 //_____________________________________________________________________________
 427 void
 428 AliMUONContourMaker::VerticalToHorizontal(const TObjArray& polygonVerticalEdges,
 429                                           TObjArray& horizontalEdges) const
 430 {
 431   /// Deduce the set of horizontal edges from the vertical edges
 432   /// Output array horizontalEdges should be empty before calling this method
 433
 434   AliCodeTimerAuto("");
 435
 436   TObjArray points; // array of AliMUONPointWithRef
 437   points.SetOwner(kTRUE);
 438
 439   SortPoints(polygonVerticalEdges,points);
 440
 441   for ( Int_t k = 0; k < (points.GetLast()+1)/2; ++k )
 442   {
 443     const AliMUONPointWithRef* p1 = static_cast<AliMUONPointWithRef*>(points.UncheckedAt(k*2));
 444     const AliMUONPointWithRef* p2 = static_cast<AliMUONPointWithRef*>(points.UncheckedAt(k*2+1));
 445
 446     const AliMUONSegment* refEdge = static_cast<const AliMUONSegment*>(polygonVerticalEdges.UncheckedAt(p1->Ref()));
 447
 448     // (p1,p2) is the horizontal edge.
 449     // refEdge is used to deduce the orientation of (p1,p2)
 450
 451     if ( AliMUONSegment::AreEqual(p1->X(),refEdge->EndX()) && AliMUONSegment::AreEqual(p1->Y(),refEdge->EndY()) )
 452 //    if ( AreEqual(p1,refEdge->End()) )
 453     {
 454       horizontalEdges.Add(new AliMUONSegment(p1->X(),p1->Y(),p2->X(),p2->Y()));
 455     }
 456     else
 457     {
 458       horizontalEdges.Add(new AliMUONSegment(p2->X(),p2->Y(),p1->X(),p1->Y()));
 459     }
 460   }
 461 }
 462