Compilation with icc
[u/mrichter/AliRoot.git] / MUON / AliMUONContourMaker.cxx
CommitLineData
0b936dc0 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
47ClassImp(AliMUONContourMaker)
48/// \endcond
49
50namespace
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//_____________________________________________________________________________
66AliMUONContourMaker::AliMUONContourMaker()
67{
cddcc1f3 68/// Default constructor
0b936dc0 69}
70
cddcc1f3 71
0b936dc0 72//_____________________________________________________________________________
73AliMUONContourMaker::~AliMUONContourMaker()
74{
cddcc1f3 75/// Destructor
0b936dc0 76}
77
78//_____________________________________________________________________________
79AliMUONContour*
80AliMUONContourMaker::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);
7da191b6 109 pol = static_cast<AliMUONPolygon*>(polygons.First());
0b936dc0 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//_____________________________________________________________________________
150AliMUONContour*
151AliMUONContourMaker::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//_____________________________________________________________________________
238void
239AliMUONContourMaker::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//_____________________________________________________________________________
261void
262AliMUONContourMaker::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//_____________________________________________________________________________
307AliMUONContour*
308AliMUONContourMaker::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//_____________________________________________________________________________
338void
339AliMUONContourMaker::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//_____________________________________________________________________________
361void
362AliMUONContourMaker::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//_____________________________________________________________________________
430void
431AliMUONContourMaker::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