Compilation with icc
[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
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