Logo ROOT  
Reference Guide
TGeoPolygon.cxx
Go to the documentation of this file.
1 // @(#)root/geom:$Id$
2 // Author: Mihaela Gheata 5/01/04
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 /** \class TGeoPolygon
13 \ingroup Geometry_classes
14 
15 An arbitrary polygon defined by vertices. The vertices
16 have to be defined CLOCKWISE in the XY plane, making either a convex
17 or concave polygon. No test for malformed polygons is performed.
18 
19 A polygon is a 2D shape defined by vertices in the XY plane. It is used by
20 TGeoXtru class for computing Contains() and Safety(). Only the pointers to
21 the actual lists of XY values are used - these are not owned by the class.
22 
23 To check if a point in XY plane is contained by a polygon, this is split
24 into an outscribed convex polygon and the remaining polygons of its subtraction
25 from the outscribed one. A point is INSIDE if it is
26 contained by the outscribed polygon but NOT by the remaining ones. Since these
27 can also be arbitrary polygons at their turn, a tree structure is formed:
28 
29 ~~~ {.cpp}
30  P = Pconvex - (Pconvex-P) where (-) means 'subtraction'
31  Pconvex-P = P1 + P2 + ... where (+) means 'union'
32 ~~~
33 
34 *Note that P1, P2, ... do not intersect each other and they are defined
35 by subsets of the list of vertices of P. They can be split in the same
36 way as P*
37 
38 Therefore, if C(P) represents the Boolean : 'does P contains a given point?',
39 then:
40 
41 ~~~ {.cpp}
42 C(P) = C(Pconvex) .and. not(C(P1) | C(P2) | ...)
43 ~~~
44 
45 For creating a polygon without TGeoXtru class, one has to call the constructor
46 TGeoPolygon(nvert) and then SetXY(Double_t *x, Double_t *y) providing the
47 arrays of X and Y vertex positions (defined clockwise) that have to 'live' longer
48 than the polygon they will describe. This complication is due to efficiency reasons.
49 At the end one has to call the FinishPolygon() method.
50 */
51 
52 #include "TGeoPolygon.h"
53 
54 #include "TObjArray.h"
55 #include "TMath.h"
56 #include "TGeoShape.h"
57 #include "TGeoManager.h"
58 #include "TVirtualGeoPainter.h"
59 
61 
62 ////////////////////////////////////////////////////////////////////////////////
63 /// Dummy constructor.
64 
66 {
67  fNvert = 0;
68  fNconvex = 0;
69  fInd = nullptr;
70  fIndc = nullptr;
71  fX = nullptr;
72  fY = nullptr;
73  fDaughters = nullptr;
76 }
77 
78 ////////////////////////////////////////////////////////////////////////////////
79 /// Default constructor.
80 
82 {
83  if (nvert<3) {
84  Fatal("Ctor", "Invalid number of vertices %i", nvert);
85  return;
86  }
87  fNvert = nvert;
88  fNconvex = 0;
89  fInd = new Int_t[nvert];
90  fIndc = nullptr;
91  fX = nullptr;
92  fY = nullptr;
93  fDaughters = nullptr;
96  SetNextIndex();
97 }
98 
99 ////////////////////////////////////////////////////////////////////////////////
100 /// Destructor
101 
103 {
104  if (fInd) delete [] fInd;
105  if (fIndc) delete [] fIndc;
106  if (fDaughters) {
107  fDaughters->Delete();
108  delete fDaughters;
109  }
110 }
111 ////////////////////////////////////////////////////////////////////////////////
112 /// Computes area of the polygon in [length^2].
113 
115 {
116  Int_t ic,i,j;
117  Double_t area = 0;
118  // Compute area of the convex part
119  for (ic=0; ic<fNvert; ic++) {
120  i = fInd[ic];
121  j = fInd[(ic+1)%fNvert];
122  area += 0.5*(fX[i]*fY[j]-fX[j]*fY[i]);
123  }
124  return TMath::Abs(area);
125 }
126 
127 ////////////////////////////////////////////////////////////////////////////////
128 /// Check if a point given by X = point[0], Y = point[1] is inside the polygon.
129 
131 {
132  Int_t i;
133  TGeoPolygon *poly;
134  for (i=0; i<fNconvex; i++)
135  if (!IsRightSided(point, fIndc[i], fIndc[(i+1)%fNconvex])) return kFALSE;
136  if (!fDaughters) return kTRUE;
138  for (i=0; i<nd; i++) {
139  poly = (TGeoPolygon*)fDaughters->UncheckedAt(i);
140  if (poly->Contains(point)) return kFALSE;
141  }
142  return kTRUE;
143 }
144 
145 ////////////////////////////////////////////////////////////////////////////////
146 /// Check polygon convexity.
147 
149 {
150  if (fNvert==3) {
151  SetConvex();
152  return;
153  }
154  Int_t j,k;
155  Double_t point[2];
156  for (Int_t i=0; i<fNvert; i++) {
157  j = (i+1)%fNvert;
158  k = (i+2)%fNvert;
159  point[0] = fX[fInd[k]];
160  point[1] = fY[fInd[k]];
161  if (!IsRightSided(point, fInd[i], fInd[j])) return;
162  }
163  SetConvex();
164 }
165 
166 ////////////////////////////////////////////////////////////////////////////////
167 /// Draw the polygon.
168 
170 {
171  if (!gGeoManager) return;
173 }
174 
175 ////////////////////////////////////////////////////////////////////////////////
176 /// Decompose polygon in a convex outscribed part and a list of daughter
177 /// polygons that have to be subtracted to get the actual one.
178 
180 {
182  // check convexity
183  ConvexCheck();
184  // find outscribed convex polygon indices
186  if (IsConvex()) {
187 // printf(" -> polygon convex -> same indices\n");
188  memcpy(fIndc, fInd, fNvert*sizeof(Int_t));
189  return;
190  }
191 // printf(" -> polygon NOT convex\n");
192  // make daughters if necessary
193  if (IsConvex()) return;
194  // ... algorithm here
195  if (!fDaughters) fDaughters = new TObjArray();
196  TGeoPolygon *poly = 0;
197  Int_t indconv = 0;
198  Int_t indnext, indback;
199  Int_t nskip;
200  while (indconv < fNconvex) {
201  indnext = (indconv+1)%fNconvex;
202  nskip = fIndc[indnext]-fIndc[indconv];
203  if (nskip<0) nskip+=fNvert;
204  if (nskip==1) {
205  indconv++;
206  continue;
207  }
208  // gap -> make polygon
209  poly = new TGeoPolygon(nskip+1);
210  poly->SetXY(fX,fY);
211  poly->SetNextIndex(fInd[fIndc[indconv]]);
212  poly->SetNextIndex(fInd[fIndc[indnext]]);
213  indback = fIndc[indnext]-1;
214  if (indback < 0) indback+=fNvert;
215  while (indback != fIndc[indconv]) {
216  poly->SetNextIndex(fInd[indback]);
217  indback--;
218  if (indback < 0) indback+=fNvert;
219  }
220  poly->FinishPolygon();
221  fDaughters->Add(poly);
222  indconv++;
223  }
224  for (indconv=0; indconv<fNconvex; indconv++) fIndc[indconv] = fInd[fIndc[indconv]];
225 }
226 
227 ////////////////////////////////////////////////////////////////////////////////
228 /// Fill list of vertices into provided arrays.
229 
231 {
232  memcpy(x, fX, fNvert*sizeof(Double_t));
233  memcpy(y, fY, fNvert*sizeof(Double_t));
234 }
235 
236 ////////////////////////////////////////////////////////////////////////////////
237 /// Fill list of vertices of the convex outscribed polygon into provided arrays.
238 
240 {
241  for (Int_t ic=0; ic<fNconvex; ic++) {
242  x[ic] = fX[fIndc[ic]];
243  y[ic] = fY[fIndc[ic]];
244  }
245 }
246 
247 ////////////////////////////////////////////////////////////////////////////////
248 /// Check if POINT is right-sided with respect to the segment defined by IND1 and IND2.
249 
250 Bool_t TGeoPolygon::IsRightSided(const Double_t *point, Int_t ind1, Int_t ind2) const
251 {
252  Double_t dot = (point[0]-fX[ind1])*(fY[ind2]-fY[ind1]) -
253  (point[1]-fY[ind1])*(fX[ind2]-fX[ind1]);
254  if (!IsClockwise()) dot = -dot;
255  if (dot<-1.E-10) return kFALSE;
256  return kTRUE;
257 }
258 
259 ////////////////////////////////////////////////////////////////////////////////
260 /// Check if a segment [0..fNvert-1] belongs to the outscribed convex pgon.
261 
263 {
264  if (i2<0) i2=(i1+1)%fNvert;
265  Double_t point[2];
266  for (Int_t i=0; i<fNvert; i++) {
267  if (i==i1 || i==i2) continue;
268  point[0] = fX[fInd[i]];
269  point[1] = fY[fInd[i]];
270  if (!IsRightSided(point, fInd[i1], fInd[i2])) return kFALSE;
271  }
272  return kTRUE;
273 }
274 
275 ////////////////////////////////////////////////////////////////////////////////
276 /// Check for illegal crossings between non-consecutive segments
277 
279 {
280  if (fNvert<4) return kFALSE;
281  Bool_t is_illegal = kFALSE;
282  Double_t x1,y1,x2,y2,x3,y3,x4,y4;
283  for (Int_t i=0; i<fNvert-2; i++) {
284  // Check segment i
285  for (Int_t j=i+2; j<fNvert; j++) {
286  // Versus segment j
287  if (i==0 && j==(fNvert-1)) continue;
288  x1 = fX[i];
289  y1 = fY[i];
290  x2 = fX[i+1];
291  y2 = fY[i+1];
292  x3 = fX[j];
293  y3 = fY[j];
294  x4 = fX[(j+1)%fNvert];
295  y4 = fY[(j+1)%fNvert];
296  if (TGeoShape::IsSegCrossing(x1,y1,x2,y2,x3,y3,x4,y4)) {
297  Error("IsIllegalCheck", "Illegal crossing of segment %d vs. segment %d", i,j);
298  is_illegal = kTRUE;
299  }
300  }
301  }
302  return is_illegal;
303 }
304 
305 ////////////////////////////////////////////////////////////////////////////////
306 /// Compute indices for the outscribed convex polygon.
307 
309 {
310  fNconvex = 0;
311  Int_t iseg = 0;
312  Int_t ivnew;
313  Bool_t conv;
314  Int_t *indconv = new Int_t[fNvert];
315  memset(indconv, 0, fNvert*sizeof(Int_t));
316  while (iseg<fNvert) {
317  if (!IsSegConvex(iseg)) {
318  if (iseg+2 > fNvert) break;
319  ivnew = (iseg+2)%fNvert;
320  conv = kFALSE;
321  // check iseg with next vertices
322  while (ivnew) {
323  if (IsSegConvex(iseg, ivnew)) {
324  conv = kTRUE;
325  break;
326  }
327  ivnew = (ivnew+1)%fNvert;
328  }
329  if (!conv) {
330 // Error("OutscribedConvex","NO convex line connection to vertex %d\n", iseg);
331  iseg++;
332  continue;
333  }
334  } else {
335  ivnew = (iseg+1)%fNvert;
336  }
337  // segment belonging to convex outscribed polygon
338  if (!fNconvex) indconv[fNconvex++] = iseg;
339  else if (indconv[fNconvex-1] != iseg) indconv[fNconvex++] = iseg;
340  if (iseg<fNvert-1) indconv[fNconvex++] = ivnew;
341  if (ivnew<iseg) break;
342  iseg = ivnew;
343  }
344  if (!fNconvex) {
345  delete [] indconv;
346  Fatal("OutscribedConvex","cannot build outscribed convex");
347  return;
348  }
349  fIndc = new Int_t[fNvert];
350  memcpy(fIndc, indconv, fNconvex*sizeof(Int_t)); // does not contain real indices yet
351  delete [] indconv;
352 }
353 
354 ////////////////////////////////////////////////////////////////////////////////
355 /// Compute minimum distance from POINT to any segment. Returns segment index.
356 
357 Double_t TGeoPolygon::Safety(const Double_t *point, Int_t &isegment) const
358 {
359  Int_t i1, i2;
360  Double_t p1[2], p2[3];
361  Double_t lsq, ssq, dx, dy, dpx, dpy, u;
362  Double_t safe=1E30;
363  Int_t isegmin=0;
364  for (i1=0; i1<fNvert; i1++) {
365  if (TGeoShape::IsSameWithinTolerance(safe,0)) {
366  isegment = isegmin;
367  return 0.;
368  }
369  i2 = (i1+1)%fNvert;
370  p1[0] = fX[i1];
371  p1[1] = fY[i1];
372  p2[0] = fX[i2];
373  p2[1] = fY[i2];
374 
375  dx = p2[0] - p1[0];
376  dy = p2[1] - p1[1];
377  dpx = point[0] - p1[0];
378  dpy = point[1] - p1[1];
379 
380  lsq = dx*dx + dy*dy;
382  ssq = dpx*dpx + dpy*dpy;
383  if (ssq < safe) {
384  safe = ssq;
385  isegmin = i1;
386  }
387  continue;
388  }
389  u = (dpx*dx + dpy*dy)/lsq;
390  if (u>1) {
391  dpx = point[0]-p2[0];
392  dpy = point[1]-p2[1];
393  } else {
394  if (u>=0) {
395  dpx -= u*dx;
396  dpy -= u*dy;
397  }
398  }
399  ssq = dpx*dpx + dpy*dpy;
400  if (ssq < safe) {
401  safe = ssq;
402  isegmin = i1;
403  }
404  }
405  isegment = isegmin;
406  safe = TMath::Sqrt(safe);
407 // printf("== segment %d: (%f, %f) - (%f, %f) safe=%f\n", isegment, fX[isegment],fY[isegment],fX[(isegment+1)%fNvert],fY[(isegment+1)%fNvert],safe);
408  return safe;
409 }
410 
411 ////////////////////////////////////////////////////////////////////////////////
412 /// Sets the next polygone index. If index<0 sets all indices consecutive
413 /// in increasing order.
414 
416 {
417  Int_t i;
418  if (index <0) {
419  for (i=0; i<fNvert; i++) fInd[i] = i;
420  return;
421  }
422  if (fNconvex >= fNvert) {
423  Error("SetNextIndex", "all indices already set");
424  return;
425  }
426  fInd[fNconvex++] = index;
427  if (fNconvex == fNvert) {
428  if (!fX || !fY) return;
429  Double_t area = 0.0;
430  for (i=0; i<fNvert; i++) area += fX[fInd[i]]*fY[fInd[(i+1)%fNvert]]-fX[fInd[(i+1)%fNvert]]*fY[fInd[i]];
431  if (area<0) TObject::SetBit(kGeoACW, kFALSE);
433  }
434 }
435 
436 ////////////////////////////////////////////////////////////////////////////////
437 /// Set X/Y array pointer for the polygon and daughters.
438 
440 {
441  Int_t i;
442  fX = x;
443  fY = y;
444  Double_t area = 0.0;
445  for (i=0; i<fNvert; i++) area += fX[fInd[i]]*fY[fInd[(i+1)%fNvert]]-fX[fInd[(i+1)%fNvert]]*fY[fInd[i]];
446  if (area<0) TObject::SetBit(kGeoACW, kFALSE);
448 
449  if (!fDaughters) return;
450  TGeoPolygon *poly;
452  for (i=0; i<nd; i++) {
453  poly = (TGeoPolygon*)fDaughters->At(i);
454  if (poly) poly->SetXY(x,y);
455  }
456 }
TGeoShape.h
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:91
TObjArray::Delete
virtual void Delete(Option_t *option="")
Remove all objects from the array AND delete all heap based objects.
Definition: TObjArray.cxx:356
TGeoPolygon::kGeoFinishPolygon
@ kGeoFinishPolygon
Definition: TGeoPolygon.h:24
TObjArray
Definition: TObjArray.h:37
x4
static const double x4[22]
Definition: RooGaussKronrodIntegrator1D.cxx:437
TGeoPolygon::Contains
Bool_t Contains(const Double_t *point) const
Check if a point given by X = point[0], Y = point[1] is inside the polygon.
Definition: TGeoPolygon.cxx:130
gGeoManager
R__EXTERN TGeoManager * gGeoManager
Definition: TGeoManager.h:602
ClassImp
#define ClassImp(name)
Definition: Rtypes.h:364
TGeoPolygon::SetNextIndex
void SetNextIndex(Int_t index=-1)
Sets the next polygone index.
Definition: TGeoPolygon.cxx:415
TGeoPolygon::GetConvexVertices
void GetConvexVertices(Double_t *x, Double_t *y) const
Fill list of vertices of the convex outscribed polygon into provided arrays.
Definition: TGeoPolygon.cxx:239
TGeoPolygon::TGeoPolygon
TGeoPolygon()
Dummy constructor.
Definition: TGeoPolygon.cxx:65
TObject::Error
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:890
TMath::Sqrt
Double_t Sqrt(Double_t x)
Definition: TMath.h:680
TGeoPolygon::ConvexCheck
void ConvexCheck()
Check polygon convexity.
Definition: TGeoPolygon.cxx:148
TGeoPolygon::IsConvex
Bool_t IsConvex() const
Definition: TGeoPolygon.h:62
TGeoPolygon::fIndc
Int_t * fIndc
Definition: TGeoPolygon.h:32
TGeoPolygon::Safety
Double_t Safety(const Double_t *point, Int_t &isegment) const
Compute minimum distance from POINT to any segment. Returns segment index.
Definition: TGeoPolygon.cxx:357
TObject::Fatal
virtual void Fatal(const char *method, const char *msgfmt,...) const
Issue fatal error message.
Definition: TObject.cxx:918
x
Double_t x[n]
Definition: legend1.C:17
TGeoPolygon::GetVertices
void GetVertices(Double_t *x, Double_t *y) const
Fill list of vertices into provided arrays.
Definition: TGeoPolygon.cxx:230
TMath::Abs
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
TObjArray::UncheckedAt
TObject * UncheckedAt(Int_t i) const
Definition: TObjArray.h:90
TObjArray::At
TObject * At(Int_t idx) const
Definition: TObjArray.h:166
bool
TObjArray::Add
void Add(TObject *obj)
Definition: TObjArray.h:74
x1
static const double x1[5]
Definition: RooGaussKronrodIntegrator1D.cxx:346
TGeoPolygon::IsIllegalCheck
Bool_t IsIllegalCheck() const
Check for illegal crossings between non-consecutive segments.
Definition: TGeoPolygon.cxx:278
TVirtualGeoPainter::DrawPolygon
virtual void DrawPolygon(const TGeoPolygon *poly)=0
Option_t
const typedef char Option_t
Definition: RtypesCore.h:66
TGeoPolygon::kGeoACW
@ kGeoACW
Definition: TGeoPolygon.h:25
x3
static const double x3[11]
Definition: RooGaussKronrodIntegrator1D.cxx:392
TGeoPolygon::IsSegConvex
Bool_t IsSegConvex(Int_t i1, Int_t i2=-1) const
Check if a segment [0..fNvert-1] belongs to the outscribed convex pgon.
Definition: TGeoPolygon.cxx:262
TGeoManager::GetGeomPainter
TVirtualGeoPainter * GetGeomPainter()
Make a default painter if none present. Returns pointer to it.
Definition: TGeoManager.cxx:2685
TGeoPolygon
Definition: TGeoPolygon.h:19
TGeoPolygon::SetXY
void SetXY(Double_t *x, Double_t *y)
Set X/Y array pointer for the polygon and daughters.
Definition: TGeoPolygon.cxx:439
TObject::SetBit
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition: TObject.cxx:696
TObjArray::GetEntriesFast
Int_t GetEntriesFast() const
Definition: TObjArray.h:64
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:92
y
Double_t y[n]
Definition: legend1.C:17
TGeoPolygon::SetConvex
void SetConvex(Bool_t flag=kTRUE)
Definition: TGeoPolygon.h:66
TGeoPolygon::IsRightSided
Bool_t IsRightSided(const Double_t *point, Int_t ind1, Int_t ind2) const
Check if POINT is right-sided with respect to the segment defined by IND1 and IND2.
Definition: TGeoPolygon.cxx:250
TGeoPolygon::fX
Double_t * fX
Definition: TGeoPolygon.h:33
TGeoPolygon::FinishPolygon
void FinishPolygon()
Decompose polygon in a convex outscribed part and a list of daughter polygons that have to be subtrac...
Definition: TGeoPolygon.cxx:179
TGeoPolygon::fInd
Int_t * fInd
Definition: TGeoPolygon.h:31
TGeoManager.h
TGeoPolygon.h
TVirtualGeoPainter.h
Double_t
double Double_t
Definition: RtypesCore.h:59
TGeoPolygon::OutscribedConvex
void OutscribedConvex()
Compute indices for the outscribed convex polygon.
Definition: TGeoPolygon.cxx:308
TObjArray.h
TGeoPolygon::fY
Double_t * fY
pointer to list of current X coordinates of vertices
Definition: TGeoPolygon.h:34
TGeoPolygon::fNconvex
Int_t fNconvex
Definition: TGeoPolygon.h:30
TGeoShape::IsSameWithinTolerance
static Bool_t IsSameWithinTolerance(Double_t a, Double_t b)
Check if two numbers differ with less than a tolerance.
Definition: TGeoShape.cxx:326
TGeoPolygon::IsClockwise
Bool_t IsClockwise() const
Definition: TGeoPolygon.h:61
x2
static const double x2[5]
Definition: RooGaussKronrodIntegrator1D.cxx:364
TGeoPolygon::fNvert
Int_t fNvert
Definition: TGeoPolygon.h:29
TGeoPolygon::fDaughters
TObjArray * fDaughters
pointer to list of current Y coordinates of vertices
Definition: TGeoPolygon.h:35
TGeoPolygon::Area
Double_t Area() const
Computes area of the polygon in [length^2].
Definition: TGeoPolygon.cxx:114
TGeoShape::IsSegCrossing
static Bool_t IsSegCrossing(Double_t x1, Double_t y1, Double_t x2, Double_t y2, Double_t x3, Double_t y3, Double_t x4, Double_t y4)
Check if segments (A,B) and (C,D) are crossing, where: A(x1,y1), B(x2,y2), C(x3,y3),...
Definition: TGeoShape.cxx:336
TMath::E
constexpr Double_t E()
Base of natural log:
Definition: TMath.h:102
TMath.h
TGeoPolygon::Draw
virtual void Draw(Option_t *option="")
Draw the polygon.
Definition: TGeoPolygon.cxx:169
TGeoPolygon::~TGeoPolygon
virtual ~TGeoPolygon()
Destructor.
Definition: TGeoPolygon.cxx:102
int