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