Logo ROOT   6.12/07
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 "TObjArray.h"
53 #include "TGeoPolygon.h"
54 #include "TMath.h"
55 #include "TGeoShape.h"
56 #include "TGeoManager.h"
57 #include "TVirtualGeoPainter.h"
58 
60 
61 ////////////////////////////////////////////////////////////////////////////////
62 /// Dummy constructor.
63 
65 {
66  fNvert = 0;
67  fNconvex = 0;
68  fInd = 0;
69  fIndc = 0;
70  fX = 0;
71  fY = 0;
72  fDaughters = 0;
75 }
76 
77 ////////////////////////////////////////////////////////////////////////////////
78 /// Default constructor.
79 
81 {
82  if (nvert<3) {
83  Fatal("Ctor", "Invalid number of vertices %i", nvert);
84  return;
85  }
86  fNvert = nvert;
87  fNconvex = 0;
88  fInd = new Int_t[nvert];
89  fIndc = 0;
90  fX = 0;
91  fY = 0;
92  fDaughters = 0;
95  SetNextIndex();
96 }
97 
98 ////////////////////////////////////////////////////////////////////////////////
99 /// Destructor
100 
102 {
103  if (fInd) delete [] fInd;
104  if (fIndc) delete [] fIndc;
105  if (fDaughters) {
106  fDaughters->Delete();
107  delete fDaughters;
108  }
109 }
110 ////////////////////////////////////////////////////////////////////////////////
111 /// Computes area of the polygon in [length^2].
112 
114 {
115  Int_t ic,i,j;
116  Double_t area = 0;
117  // Compute area of the convex part
118  for (ic=0; ic<fNvert; ic++) {
119  i = fInd[ic];
120  j = fInd[(ic+1)%fNvert];
121  area += 0.5*(fX[i]*fY[j]-fX[j]*fY[i]);
122  }
123  return TMath::Abs(area);
124 }
125 
126 ////////////////////////////////////////////////////////////////////////////////
127 /// Check if a point given by X = point[0], Y = point[1] is inside the polygon.
128 
130 {
131  Int_t i;
132  TGeoPolygon *poly;
133  for (i=0; i<fNconvex; i++)
134  if (!IsRightSided(point, fIndc[i], fIndc[(i+1)%fNconvex])) return kFALSE;
135  if (!fDaughters) return kTRUE;
137  for (i=0; i<nd; i++) {
138  poly = (TGeoPolygon*)fDaughters->UncheckedAt(i);
139  if (poly->Contains(point)) return kFALSE;
140  }
141  return kTRUE;
142 }
143 
144 ////////////////////////////////////////////////////////////////////////////////
145 /// Check polygon convexity.
146 
148 {
149  if (fNvert==3) {
150  SetConvex();
151  return;
152  }
153  Int_t j,k;
154  Double_t point[2];
155  for (Int_t i=0; i<fNvert; i++) {
156  j = (i+1)%fNvert;
157  k = (i+2)%fNvert;
158  point[0] = fX[fInd[k]];
159  point[1] = fY[fInd[k]];
160  if (!IsRightSided(point, fInd[i], fInd[j])) return;
161  }
162  SetConvex();
163 }
164 
165 ////////////////////////////////////////////////////////////////////////////////
166 /// Draw the polygon.
167 
169 {
170  if (!gGeoManager) return;
172 }
173 
174 ////////////////////////////////////////////////////////////////////////////////
175 /// Decompose polygon in a convex outscribed part and a list of daughter
176 /// polygons that have to be subtracted to get the actual one.
177 
179 {
181  // check convexity
182  ConvexCheck();
183  // find outscribed convex polygon indices
185  if (IsConvex()) {
186 // printf(" -> polygon convex -> same indices\n");
187  memcpy(fIndc, fInd, fNvert*sizeof(Int_t));
188  return;
189  }
190 // printf(" -> polygon NOT convex\n");
191  // make daughters if necessary
192  if (IsConvex()) return;
193  // ... algorithm here
194  if (!fDaughters) fDaughters = new TObjArray();
195  TGeoPolygon *poly = 0;
196  Int_t indconv = 0;
197  Int_t indnext, indback;
198  Int_t nskip;
199  while (indconv < fNconvex) {
200  indnext = (indconv+1)%fNconvex;
201  nskip = fIndc[indnext]-fIndc[indconv];
202  if (nskip<0) nskip+=fNvert;
203  if (nskip==1) {
204  indconv++;
205  continue;
206  }
207  // gap -> make polygon
208  poly = new TGeoPolygon(nskip+1);
209  poly->SetXY(fX,fY);
210  poly->SetNextIndex(fInd[fIndc[indconv]]);
211  poly->SetNextIndex(fInd[fIndc[indnext]]);
212  indback = fIndc[indnext]-1;
213  if (indback < 0) indback+=fNvert;
214  while (indback != fIndc[indconv]) {
215  poly->SetNextIndex(fInd[indback]);
216  indback--;
217  if (indback < 0) indback+=fNvert;
218  }
219  poly->FinishPolygon();
220  fDaughters->Add(poly);
221  indconv++;
222  }
223  for (indconv=0; indconv<fNconvex; indconv++) fIndc[indconv] = fInd[fIndc[indconv]];
224 }
225 
226 ////////////////////////////////////////////////////////////////////////////////
227 /// Fill list of vertices into provided arrays.
228 
230 {
231  memcpy(x, fX, fNvert*sizeof(Double_t));
232  memcpy(y, fY, fNvert*sizeof(Double_t));
233 }
234 
235 ////////////////////////////////////////////////////////////////////////////////
236 /// Fill list of vertices of the convex outscribed polygon into provided arrays.
237 
239 {
240  for (Int_t ic=0; ic<fNconvex; ic++) {
241  x[ic] = fX[fIndc[ic]];
242  y[ic] = fY[fIndc[ic]];
243  }
244 }
245 
246 ////////////////////////////////////////////////////////////////////////////////
247 /// Check if POINT is right-sided with respect to the segment defined by IND1 and IND2.
248 
249 Bool_t TGeoPolygon::IsRightSided(const Double_t *point, Int_t ind1, Int_t ind2) const
250 {
251  Double_t dot = (point[0]-fX[ind1])*(fY[ind2]-fY[ind1]) -
252  (point[1]-fY[ind1])*(fX[ind2]-fX[ind1]);
253  if (!IsClockwise()) dot = -dot;
254  if (dot<-1.E-10) return kFALSE;
255  return kTRUE;
256 }
257 
258 ////////////////////////////////////////////////////////////////////////////////
259 /// Check if a segment [0..fNvert-1] belongs to the outscribed convex pgon.
260 
262 {
263  if (i2<0) i2=(i1+1)%fNvert;
264  Double_t point[2];
265  for (Int_t i=0; i<fNvert; i++) {
266  if (i==i1 || i==i2) continue;
267  point[0] = fX[fInd[i]];
268  point[1] = fY[fInd[i]];
269  if (!IsRightSided(point, fInd[i1], fInd[i2])) return kFALSE;
270  }
271  return kTRUE;
272 }
273 
274 ////////////////////////////////////////////////////////////////////////////////
275 /// Check for illegal crossings between non-consecutive segments
276 
278 {
279  if (fNvert<4) return kFALSE;
280  Bool_t is_illegal = kFALSE;
281  Double_t x1,y1,x2,y2,x3,y3,x4,y4;
282  for (Int_t i=0; i<fNvert-2; i++) {
283  // Check segment i
284  for (Int_t j=i+2; j<fNvert; j++) {
285  // Versus segment j
286  if (i==0 && j==(fNvert-1)) continue;
287  x1 = fX[i];
288  y1 = fY[i];
289  x2 = fX[i+1];
290  y2 = fY[i+1];
291  x3 = fX[j];
292  y3 = fY[j];
293  x4 = fX[(j+1)%fNvert];
294  y4 = fY[(j+1)%fNvert];
295  if (TGeoShape::IsSegCrossing(x1,y1,x2,y2,x3,y3,x4,y4)) {
296  Error("IsIllegalCheck", "Illegal crossing of segment %d vs. segment %d", i,j);
297  is_illegal = kTRUE;
298  }
299  }
300  }
301  return is_illegal;
302 }
303 
304 ////////////////////////////////////////////////////////////////////////////////
305 /// Compute indices for the outscribed convex polygon.
306 
308 {
309  fNconvex = 0;
310  Int_t iseg = 0;
311  Int_t ivnew;
312  Bool_t conv;
313  Int_t *indconv = new Int_t[fNvert];
314  memset(indconv, 0, fNvert*sizeof(Int_t));
315  while (iseg<fNvert) {
316  if (!IsSegConvex(iseg)) {
317  if (iseg+2 > fNvert) break;
318  ivnew = (iseg+2)%fNvert;
319  conv = kFALSE;
320  // check iseg with next vertices
321  while (ivnew) {
322  if (IsSegConvex(iseg, ivnew)) {
323  conv = kTRUE;
324  break;
325  }
326  ivnew = (ivnew+1)%fNvert;
327  }
328  if (!conv) {
329 // Error("OutscribedConvex","NO convex line connection to vertex %d\n", iseg);
330  iseg++;
331  continue;
332  }
333  } else {
334  ivnew = (iseg+1)%fNvert;
335  }
336  // segment belonging to convex outscribed polygon
337  if (!fNconvex) indconv[fNconvex++] = iseg;
338  else if (indconv[fNconvex-1] != iseg) indconv[fNconvex++] = iseg;
339  if (iseg<fNvert-1) indconv[fNconvex++] = ivnew;
340  if (ivnew<iseg) break;
341  iseg = ivnew;
342  }
343  if (!fNconvex) {
344  delete [] indconv;
345  Fatal("OutscribedConvex","cannot build outscribed convex");
346  return;
347  }
348  fIndc = new Int_t[fNvert];
349  memcpy(fIndc, indconv, fNconvex*sizeof(Int_t)); // does not contain real indices yet
350  delete [] indconv;
351 }
352 
353 ////////////////////////////////////////////////////////////////////////////////
354 /// Compute minimum distance from POINT to any segment. Returns segment index.
355 
356 Double_t TGeoPolygon::Safety(const Double_t *point, Int_t &isegment) const
357 {
358  Int_t i1, i2;
359  Double_t p1[2], p2[3];
360  Double_t lsq, ssq, dx, dy, dpx, dpy, u;
361  Double_t safe=1E30;
362  Int_t isegmin=0;
363  for (i1=0; i1<fNvert; i1++) {
364  if (TGeoShape::IsSameWithinTolerance(safe,0)) {
365  isegment = isegmin;
366  return 0.;
367  }
368  i2 = (i1+1)%fNvert;
369  p1[0] = fX[i1];
370  p1[1] = fY[i1];
371  p2[0] = fX[i2];
372  p2[1] = fY[i2];
373 
374  dx = p2[0] - p1[0];
375  dy = p2[1] - p1[1];
376  dpx = point[0] - p1[0];
377  dpy = point[1] - p1[1];
378 
379  lsq = dx*dx + dy*dy;
381  ssq = dpx*dpx + dpy*dpy;
382  if (ssq < safe) {
383  safe = ssq;
384  isegmin = i1;
385  }
386  continue;
387  }
388  u = (dpx*dx + dpy*dy)/lsq;
389  if (u>1) {
390  dpx = point[0]-p2[0];
391  dpy = point[1]-p2[1];
392  } else {
393  if (u>=0) {
394  dpx -= u*dx;
395  dpy -= u*dy;
396  }
397  }
398  ssq = dpx*dpx + dpy*dpy;
399  if (ssq < safe) {
400  safe = ssq;
401  isegmin = i1;
402  }
403  }
404  isegment = isegmin;
405  safe = TMath::Sqrt(safe);
406 // printf("== segment %d: (%f, %f) - (%f, %f) safe=%f\n", isegment, fX[isegment],fY[isegment],fX[(isegment+1)%fNvert],fY[(isegment+1)%fNvert],safe);
407  return safe;
408 }
409 
410 ////////////////////////////////////////////////////////////////////////////////
411 /// Sets the next polygone index. If index<0 sets all indices consecutive
412 /// in increasing order.
413 
415 {
416  Int_t i;
417  if (index <0) {
418  for (i=0; i<fNvert; i++) fInd[i] = i;
419  return;
420  }
421  if (fNconvex >= fNvert) {
422  Error("SetNextIndex", "all indices already set");
423  return;
424  }
425  fInd[fNconvex++] = index;
426  if (fNconvex == fNvert) {
427  if (!fX || !fY) return;
428  Double_t area = 0.0;
429  for (i=0; i<fNvert; i++) area += fX[fInd[i]]*fY[fInd[(i+1)%fNvert]]-fX[fInd[(i+1)%fNvert]]*fY[fInd[i]];
430  if (area<0) TObject::SetBit(kGeoACW, kFALSE);
432  }
433 }
434 
435 ////////////////////////////////////////////////////////////////////////////////
436 /// Set X/Y array pointer for the polygon and daughters.
437 
439 {
440  Int_t i;
441  fX = x;
442  fY = y;
443  Double_t area = 0.0;
444  for (i=0; i<fNvert; i++) area += fX[fInd[i]]*fY[fInd[(i+1)%fNvert]]-fX[fInd[(i+1)%fNvert]]*fY[fInd[i]];
445  if (area<0) TObject::SetBit(kGeoACW, kFALSE);
447 
448  if (!fDaughters) return;
449  TGeoPolygon *poly;
451  for (i=0; i<nd; i++) {
452  poly = (TGeoPolygon*)fDaughters->At(i);
453  if (poly) poly->SetXY(x,y);
454  }
455 }
virtual void Draw(Option_t *option="")
Draw the polygon.
An array of TObjects.
Definition: TObjArray.h:37
void SetNextIndex(Int_t index=-1)
Sets the next polygone index.
const char Option_t
Definition: RtypesCore.h:62
virtual void Delete(Option_t *option="")
Remove all objects from the array AND delete all heap based objects.
Definition: TObjArray.cxx:355
void ConvexCheck()
Check polygon convexity.
Bool_t Contains(const Double_t *point) const
Check if a point given by X = point[0], Y = point[1] is inside the polygon.
Bool_t IsConvex() const
Definition: TGeoPolygon.h:59
Int_t * fIndc
Definition: TGeoPolygon.h:32
Double_t Safety(const Double_t *point, Int_t &isegment) const
Compute minimum distance from POINT to any segment. Returns segment index.
An arbitrary polygon defined by vertices.
Definition: TGeoPolygon.h:19
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
void GetConvexVertices(Double_t *x, Double_t *y) const
Fill list of vertices of the convex outscribed polygon into provided arrays.
void GetVertices(Double_t *x, Double_t *y) const
Fill list of vertices into provided arrays.
TGeoPolygon()
Dummy constructor.
Definition: TGeoPolygon.cxx:64
TObject * At(Int_t idx) const
Definition: TObjArray.h:165
Short_t Abs(Short_t d)
Definition: TMathBase.h:108
static Bool_t IsSameWithinTolerance(Double_t a, Double_t b)
Check if two numbers differ with less than a tolerance.
Definition: TGeoShape.cxx:326
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition: TObject.cxx:694
static const double x2[5]
Double_t x[n]
Definition: legend1.C:17
Bool_t IsIllegalCheck() const
Check for illegal crossings between non-consecutive segments.
static double p2(double t, double a, double b, double c)
static const double x4[22]
Bool_t IsSegConvex(Int_t i1, Int_t i2=-1) const
Check if a segment [0..fNvert-1] belongs to the outscribed convex pgon.
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), D(x4,y4)
Definition: TGeoShape.cxx:336
void SetXY(Double_t *x, Double_t *y)
Set X/Y array pointer for the polygon and daughters.
Int_t GetEntriesFast() const
Definition: TObjArray.h:64
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:880
static double p1(double t, double a, double b)
Double_t * fX
Definition: TGeoPolygon.h:33
constexpr Double_t E()
Definition: TMath.h:74
const Bool_t kFALSE
Definition: RtypesCore.h:88
void SetConvex(Bool_t flag=kTRUE)
Definition: TGeoPolygon.h:63
void OutscribedConvex()
Compute indices for the outscribed convex polygon.
TObject * UncheckedAt(Int_t i) const
Definition: TObjArray.h:89
static const double x1[5]
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.
#define ClassImp(name)
Definition: Rtypes.h:359
R__EXTERN TGeoManager * gGeoManager
Definition: TGeoManager.h:559
double Double_t
Definition: RtypesCore.h:55
Double_t * fY
pointer to list of current X coordinates of vertices
Definition: TGeoPolygon.h:34
void FinishPolygon()
Decompose polygon in a convex outscribed part and a list of daughter polygons that have to be subtrac...
Int_t * fInd
Definition: TGeoPolygon.h:31
Bool_t IsClockwise() const
Definition: TGeoPolygon.h:58
Double_t y[n]
Definition: legend1.C:17
virtual void DrawPolygon(const TGeoPolygon *poly)=0
TObjArray * fDaughters
pointer to list of current Y coordinates of vertices
Definition: TGeoPolygon.h:35
Int_t fNconvex
Definition: TGeoPolygon.h:30
void Add(TObject *obj)
Definition: TObjArray.h:73
Double_t Sqrt(Double_t x)
Definition: TMath.h:590
Int_t fNvert
Definition: TGeoPolygon.h:29
virtual void Fatal(const char *method, const char *msgfmt,...) const
Issue fatal error message.
Definition: TObject.cxx:908
Double_t Area() const
Computes area of the polygon in [length^2].
virtual ~TGeoPolygon()
Destructor.
const Bool_t kTRUE
Definition: RtypesCore.h:87
TVirtualGeoPainter * GetGeomPainter()
Make a default painter if none present. Returns pointer to it.
static const double x3[11]