Adding classes forgotten in previous commit
[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 }
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