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 /// Default constructor
  69 }
  70
  71
  72 //_____________________________________________________________________________
  73 AliMUONContourMaker::~AliMUONContourMaker()
  74 {
  75 /// Destructor
  76 }
  77
  78 //_____________________________________________________________________________
  79 AliMUONContour*
  80 AliMUONContourMaker::CreateContour(const TObjArray& polygons, const char* name) const
  81 {
  82   /// Create the contour of the polygon array
  83   /// and get back the intermediate verticals and horizontal segments
  84   /// both arrays are arrays of AliMUONSegment objects.
  85
  86   AliCodeTimerAuto("");
  87
  88   if ( polygons.IsEmpty() ) return 0x0; // protection against user error...
  89
  90   // Sanity check : insure that all polygons are oriented counter-clockwise
  91   TIter next(&polygons);
  92   AliMUONPolygon* pol;
  93   while ( ( pol = static_cast<AliMUONPolygon*>(next()) ) )
  94   {
  95     if ( !pol->IsCounterClockwiseOriented() )
  96     {
  97       AliError(Form("Got a clockwise oriented polygon in CreateContour(%s). That's not OK !",name));
  98       StdoutToAliError(polygons.Print());
  99       return 0x0;
 100     }
 101   }
 102
 103   AliMUONContour* contour(0x0);
 104
 105   if ( polygons.GetLast() == 0 )
 106   {
 107     AliCodeTimerAuto("Trivial case");
 108     contour = new AliMUONContour(name);
 109     pol = static_cast<AliMUONPolygon*>(polygons.First());
 110     contour->Add(*pol);
 111     contour->AssertOrientation();
 112     return contour;
 113   }
 114
 115   TObjArray polygonVerticalEdges; // arrray of AliMUONSegment objects
 116   polygonVerticalEdges.SetOwner(kTRUE);
 117   // get vertical edges of input polygons
 118   GetVerticalEdges(polygons,polygonVerticalEdges);
 119
 120   // sort them in ascending x order
 121   // if same x, insure that left edges are before right edges
 122   // within same x, order by increasing bottommost y (see AliMUONSegment::Compare method)
 123   polygonVerticalEdges.Sort();
 124
 125   if ( polygonVerticalEdges.GetLast()+1 < 2 )
 126   {
 127     AliError(Form("Got too few edges here for createContour %s",name));
 128     polygons.Print();
 129     TObject* o(0x0);
 130     o->Print();
 131   }
 132
 133   // Find the vertical edges of the merged contour. This is the meat of the algorithm...
 134   TObjArray contourVerticalEdges;
 135   contourVerticalEdges.SetOwner(kTRUE);
 136   Sweep(polygonVerticalEdges,contourVerticalEdges);
 137
 138   TObjArray horizontals;
 139   horizontals.SetOwner(kTRUE);
 140   VerticalToHorizontal(contourVerticalEdges,horizontals);
 141
 142   contour = FinalizeContour(contourVerticalEdges,horizontals);
 143
 144   if ( contour && name ) contour->SetName(name);
 145
 146   return contour;
 147 }
 148
 149 //_____________________________________________________________________________
 150 AliMUONContour*
 151 AliMUONContourMaker::FinalizeContour(const TObjArray& verticals,
 152                                      const TObjArray& horizontals) const
 153 {
 154   /// For a list of vertical and horizontal edges, we build the final
 155   /// contour object.
 156
 157   AliCodeTimerAuto("");
 158
 159   TObjArray all; // array of AliMUONSegment
 160   TObjArray inorder; // array of AliMUONSegment
 161
 162   all.SetOwner(kFALSE);
 163   inorder.SetOwner(kFALSE);
 164
 165   for ( Int_t i = 0; i <= verticals.GetLast(); ++i )
 166   {
 167     all.Add(verticals.UncheckedAt(i));
 168     all.Add(horizontals.UncheckedAt(i));
 169   }
 170
 171   Int_t i(0);
 172
 173   AliMUONContour* contour = new AliMUONContour;
 174
 175   int total(0);
 176
 177   while ( !all.IsEmpty() )
 178   {
 179     total++;
 180
 181     if ( total > 1000 )
 182     {
 183       AliError("Total 1000 reached !!!!");
 184       return 0x0;
 185     }
 186
 187     AliMUONSegment* si = static_cast<AliMUONSegment*>(all.UncheckedAt(i));
 188     inorder.Add(si);
 189     const AliMUONSegment* all0 = static_cast<const AliMUONSegment*>(all.First());
 190     if ( i != 0 && AliMUONSegment::AreEqual(si->EndX(),all0->StartX()) && AliMUONSegment::AreEqual(si->EndY(),all0->StartY()) )
 191     {
 192       Int_t n(-1);
 193
 194       AliMUONPolygon polygon(inorder.GetLast()+2);
 195
 196       // we got a cycle. Add it to the contour
 197       for ( Int_t j = 0; j <= inorder.GetLast(); ++j )
 198       {
 199         AliMUONSegment* s = static_cast<AliMUONSegment*>(inorder.UncheckedAt(j));
 200         polygon.SetVertex(++n,s->StartX(),s->StartY());
 201         all.Remove(s);
 202       }
 203
 204       all.Compress();
 205
 206       polygon.Close();
 207
 208       contour->Add(polygon);
 209
 210       if ( ! all.IsEmpty() )
 211       {
 212         i = 0;
 213         inorder.Clear();
 214       }
 215       continue;
 216     }
 217
 218     for ( Int_t j = 0; j <= all.GetLast(); ++j)
 219     {
 220       if ( j != i )
 221       {
 222         const AliMUONSegment* sj = static_cast<const AliMUONSegment*>(all.UncheckedAt(j));
 223         if ( AliMUONSegment::AreEqual(si->EndX(),sj->StartX()) && AliMUONSegment::AreEqual(si->EndY(),sj->StartY()))
 224         {
 225           i = j;
 226           break;
 227         }
 228       }
 229     }
 230   }
 231
 232   contour->AssertOrientation(kTRUE);
 233   return contour;
 234 }
 235
 236
 237 //_____________________________________________________________________________
 238 void
 239 AliMUONContourMaker::GetVerticalEdges(const TObjArray& polygons, TObjArray& polygonVerticalEdges) const
 240 {
 241   /// From an array of polygons, extract the list of vertical edges.
 242   /// Output array polygonVerticalEdges should be empty before calling.
 243
 244   AliCodeTimerAuto("");
 245
 246   for ( Int_t i = 0; i <= polygons.GetLast(); ++i )
 247   {
 248     const AliMUONPolygon* g = static_cast<const AliMUONPolygon*>(polygons.UncheckedAt(i));
 249     for ( Int_t j = 0; j < g->NumberOfVertices()-1; ++j )
 250     {
 251       if ( AliMUONSegment::AreEqual(g->X(j),g->X(j+1)) ) // segment is vertical
 252       {
 253         polygonVerticalEdges.Add(new AliMUONSegment(g->X(j),g->Y(j),g->X(j+1),g->Y(j+1)));
 254       }
 255     }
 256   }
 257 }
 258
 259
 260 //_____________________________________________________________________________
 261 void
 262 AliMUONContourMaker::GetYPositions(const TObjArray& polygonVerticalEdges,
 263                                    TArrayD& yPositions) const
 264 {
 265   /// Fill the array yPositions with the different y positions found in
 266   /// polygonVerticalEdges
 267
 268   AliCodeTimerAuto("");
 269
 270   Double_t* y = new Double_t[polygonVerticalEdges.GetSize()*2];
 271   Int_t n(0);
 272
 273   for ( Int_t i = 0; i < polygonVerticalEdges.GetLast(); ++i )
 274   {
 275     AliMUONSegment* s = static_cast<AliMUONSegment*>(polygonVerticalEdges.UncheckedAt(i));
 276     y[n] = s->StartY();
 277     y[n+1] = s->EndY();
 278     n += 2;
 279   }
 280   Int_t* ix = new Int_t[n+1];
 281
 282   TMath::Sort(n,y,ix,kFALSE);
 283
 284   yPositions.Set(n+1);
 285
 286   Int_t u(0);
 287   Double_t x(FLT_MAX);
 288
 289   for ( Int_t i = 0; i < n; ++i )
 290   {
 291     if ( y[ix[i]] != x )
 292     {
 293       yPositions[u] = y[ix[i]];
 294       x = y[ix[i]];
 295       ++u;
 296     }
 297   }
 298
 299   yPositions.Set(u);
 300
 301   delete[] ix;
 302   delete[] y;
 303
 304 }
 305
 306 //_____________________________________________________________________________
 307 AliMUONContour*
 308 AliMUONContourMaker::MergeContour(const TObjArray& contours, const char* name) const
 309 {
 310   /// Merge all the polygons of all contours into a single contour
 311
 312   AliCodeTimerAuto("");
 313
 314   TObjArray polygons;
 315   polygons.SetOwner(kTRUE);
 316
 317   TIter next(&contours);
 318   AliMUONContour* contour;
 319   while ( ( contour = static_cast<AliMUONContour*>(next()) ) )
 320   {
 321     const TObjArray* contourPolygons = contour->Polygons();
 322     TIter nextPol(contourPolygons);
 323     AliMUONPolygon* pol;
 324     while ( ( pol = static_cast<AliMUONPolygon*>(nextPol()) ) )
 325     {
 326       polygons.Add(new AliMUONPolygon(*pol));
 327     }
 328   }
 329
 330   if ( polygons.IsEmpty() ) return 0x0;
 331
 332   contour = CreateContour(polygons,name);
 333
 334   return contour;
 335 }
 336
 337 //_____________________________________________________________________________
 338 void
 339 AliMUONContourMaker::SortPoints(const TObjArray& polygonVerticalEdges,
 340                                 TObjArray& sortedPoints) const
 341 {
 342   /// Sort the point of the vertical edges in ascending order, first on ordinate,
 343   /// then on abcissa, and put them in output vector sortedPoints.
 344   /// Output array sortedPoints should be empty before calling this method.
 345
 346   AliCodeTimerAuto("");
 347
 348   for ( Int_t i = 0; i <= polygonVerticalEdges.GetLast(); ++i )
 349   {
 350     const AliMUONSegment* e = static_cast<const AliMUONSegment*>(polygonVerticalEdges.UncheckedAt(i));
 351     sortedPoints.Add(new AliMUONPointWithRef(e->StartX(),e->StartY(),i));
 352     sortedPoints.Add(new AliMUONPointWithRef(e->EndX(),e->EndY(),i));
 353     // note that we keep track of the original edge, which is used
 354     // later on to deduce orientation of horizontal edges.
 355   }
 356
 357   sortedPoints.Sort();
 358 }
 359
 360 //_____________________________________________________________________________
 361 void
 362 AliMUONContourMaker::Sweep(const TObjArray& polygonVerticalEdges,
 363                            TObjArray& contourVerticalEdges) const
 364 {
 365   /// This is the meat of the algorithm of the contour merging...
 366
 367   AliCodeTimerAuto("");
 368
 369   TArrayD yPositions;
 370   GetYPositions(polygonVerticalEdges,yPositions);
 371
 372   AliMUONSegmentTree segmentTree(yPositions);
 373
 374   for ( Int_t i = 0; i <= polygonVerticalEdges.GetLast(); ++i )
 375   {
 376     const AliMUONSegment* edge = static_cast<const AliMUONSegment*>(polygonVerticalEdges.UncheckedAt(i));
 377
 378     assert(edge!=0x0);
 379
 380     if ( edge->IsLeftEdge() )
 381     {
 382       segmentTree.Contribution(edge->Bottom(),edge->Top());
 383       segmentTree.InsertInterval(edge->Bottom(),edge->Top());
 384     }
 385     else
 386     {
 387       segmentTree.DeleteInterval(edge->Bottom(),edge->Top());
 388       segmentTree.Contribution(edge->Bottom(),edge->Top());
 389     }
 390
 391     AliMUONSegment e1(*edge);
 392
 393     if ( i < polygonVerticalEdges.GetLast() )
 394     {
 395       const AliMUONSegment* next = static_cast<const AliMUONSegment*>(polygonVerticalEdges.UncheckedAt(i+1));
 396       e1 = *next;
 397     }
 398
 399     if ( ( edge->IsLeftEdge() != e1.IsLeftEdge() ) ||
 400         ( !AliMUONSegment::AreEqual(edge->StartX(),e1.StartX() ) ) ||
 401         ( i == polygonVerticalEdges.GetLast() ) )
 402     {
 403       const TObjArray& stack = segmentTree.Stack();
 404
 405       double x = edge->StartX();
 406
 407       for ( Int_t j = 0; j <= stack.GetLast(); ++j )
 408       {
 409         AliMUONSegment* sj = static_cast<AliMUONSegment*>(stack.UncheckedAt(j));
 410         AliMUONSegment* s = new AliMUONSegment(x,sj->StartY(),x,sj->EndY());
 411
 412         if  (s->IsAPoint())
 413         {
 414           delete s;
 415           continue;
 416         }
 417
 418         if ( edge->IsLeftEdge() != s->IsLeftEdge() )
 419         {
 420           s->Set(x,sj->EndY(),x,sj->StartY());
 421         }
 422         contourVerticalEdges.Add(s);
 423       }
 424       segmentTree.ResetStack();
 425     }
 426   }
 427 }
 428
 429 //_____________________________________________________________________________
 430 void
 431 AliMUONContourMaker::VerticalToHorizontal(const TObjArray& polygonVerticalEdges,
 432                                           TObjArray& horizontalEdges) const
 433 {
 434   /// Deduce the set of horizontal edges from the vertical edges
 435   /// Output array horizontalEdges should be empty before calling this method
 436
 437   AliCodeTimerAuto("");
 438
 439   TObjArray points; // array of AliMUONPointWithRef
 440   points.SetOwner(kTRUE);
 441
 442   SortPoints(polygonVerticalEdges,points);
 443
 444   for ( Int_t k = 0; k < (points.GetLast()+1)/2; ++k )
 445   {
 446     const AliMUONPointWithRef* p1 = static_cast<AliMUONPointWithRef*>(points.UncheckedAt(k*2));
 447     const AliMUONPointWithRef* p2 = static_cast<AliMUONPointWithRef*>(points.UncheckedAt(k*2+1));
 448
 449     const AliMUONSegment* refEdge = static_cast<const AliMUONSegment*>(polygonVerticalEdges.UncheckedAt(p1->Ref()));
 450
 451     // (p1,p2) is the horizontal edge.
 452     // refEdge is used to deduce the orientation of (p1,p2)
 453
 454     if ( AliMUONSegment::AreEqual(p1->X(),refEdge->EndX()) && AliMUONSegment::AreEqual(p1->Y(),refEdge->EndY()) )
 455 //    if ( AreEqual(p1,refEdge->End()) )
 456     {
 457       horizontalEdges.Add(new AliMUONSegment(p1->X(),p1->Y(),p2->X(),p2->Y()));
 458     }
 459     else
 460     {
 461       horizontalEdges.Add(new AliMUONSegment(p2->X(),p2->Y(),p1->X(),p1->Y()));
 462     }
 463   }
 464 }
 465