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