 ROOT   Reference Guide
Searching...
No Matches
TGeoPolygon.cxx
Go to the documentation of this file.
1// @(#)root/geom:$Id$
2// Author: Mihaela Gheata 5/01/04
3
4/*************************************************************************
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
15An arbitrary polygon defined by vertices. The vertices
16have to be defined CLOCKWISE in the XY plane, making either a convex
17or concave polygon. No test for malformed polygons is performed.
18
19A polygon is a 2D shape defined by vertices in the XY plane. It is used by
20TGeoXtru class for computing Contains() and Safety(). Only the pointers to
21the actual lists of XY values are used - these are not owned by the class.
22
23To check if a point in XY plane is contained by a polygon, this is split
24into an outscribed convex polygon and the remaining polygons of its subtraction
25from the outscribed one. A point is INSIDE if it is
26contained by the outscribed polygon but NOT by the remaining ones. Since these
27can 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
35by subsets of the list of vertices of P. They can be split in the same
36way as P*
37
38Therefore, if C(P) represents the Boolean : 'does P contains a given point?',
39then:
40
41~~~ {.cpp}
42C(P) = C(Pconvex) .and. not(C(P1) | C(P2) | ...)
43~~~
44
45For creating a polygon without TGeoXtru class, one has to call the constructor
46TGeoPolygon(nvert) and then SetXY(Double_t *x, Double_t *y) providing the
47arrays of X and Y vertex positions (defined clockwise) that have to 'live' longer
48than the polygon they will describe. This complication is due to efficiency reasons.
49At 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;
97}
98
99////////////////////////////////////////////////////////////////////////////////
100/// Destructor
101
103{
104 if (fInd) delete [] fInd;
105 if (fIndc) delete [] fIndc;
106 if (fDaughters) {
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, Y = point 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++) {
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;
156 for (Int_t i=0; i<fNvert; i++) {
157 j = (i+1)%fNvert;
158 k = (i+2)%fNvert;
159 point = fX[fInd[k]];
160 point = 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();
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
250Bool_t TGeoPolygon::IsRightSided(const Double_t *point, Int_t ind1, Int_t ind2) const
251{
252 Double_t dot = (point-fX[ind1])*(fY[ind2]-fY[ind1]) -
253 (point-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;
266 for (Int_t i=0; i<fNvert; i++) {
267 if (i==i1 || i==i2) continue;
268 point = fX[fInd[i]];
269 point = 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
357Double_t TGeoPolygon::Safety(const Double_t *point, Int_t &isegment) const
358{
359 Int_t i1, i2;
360 Double_t p1, p2;
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++) {
366 isegment = isegmin;
367 return 0.;
368 }
369 i2 = (i1+1)%fNvert;
370 p1 = fX[i1];
371 p1 = fY[i1];
372 p2 = fX[i2];
373 p2 = fY[i2];
374
375 dx = p2 - p1;
376 dy = p2 - p1;
377 dpx = point - p1;
378 dpy = point - p1;
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-p2;
392 dpy = point-p2;
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}
static const double x2
static const double x4
static const double x1
static const double x3
const Bool_t kFALSE
Definition RtypesCore.h:92
double Double_t
Definition RtypesCore.h:59
const Bool_t kTRUE
Definition RtypesCore.h:91
const char Option_t
Definition RtypesCore.h:66
#define ClassImp(name)
Definition Rtypes.h:364
R__EXTERN TGeoManager * gGeoManager
TVirtualGeoPainter * GetGeomPainter()
Make a default painter if none present. Returns pointer to it.
An arbitrary polygon defined by vertices.
Definition TGeoPolygon.h:20
Bool_t IsConvex() const
Definition TGeoPolygon.h:62
Bool_t IsSegConvex(Int_t i1, Int_t i2=-1) const
Check if a segment [0..fNvert-1] belongs to the outscribed convex pgon.
TGeoPolygon()
Dummy constructor.
Double_t Safety(const Double_t *point, Int_t &isegment) const
Compute minimum distance from POINT to any segment. Returns segment index.
void SetXY(Double_t *x, Double_t *y)
Set X/Y array pointer for the polygon and daughters.
Double_t Area() const
Computes area of the polygon in [length^2].
Bool_t Contains(const Double_t *point) const
Check if a point given by X = point, Y = point is inside the polygon.
Bool_t IsIllegalCheck() const
Check for illegal crossings between non-consecutive segments.
Int_t fNconvex
Definition TGeoPolygon.h:30
void GetVertices(Double_t *x, Double_t *y) const
Fill list of vertices into provided arrays.
void SetNextIndex(Int_t index=-1)
Sets the next polygone index.
Int_t * fInd
Definition TGeoPolygon.h:31
virtual void Draw(Option_t *option="")
Draw the polygon.
virtual ~TGeoPolygon()
Destructor.
Double_t * fY
pointer to list of current X coordinates of vertices
Definition TGeoPolygon.h:34
Bool_t IsClockwise() const
Definition TGeoPolygon.h:61
void OutscribedConvex()
Compute indices for the outscribed convex polygon.
Double_t * fX
Definition TGeoPolygon.h:33
Int_t * fIndc
Definition TGeoPolygon.h:32
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.
void ConvexCheck()
Check polygon convexity.
TObjArray * fDaughters
pointer to list of current Y coordinates of vertices
Definition TGeoPolygon.h:35
void GetConvexVertices(Double_t *x, Double_t *y) const
Fill list of vertices of the convex outscribed polygon into provided arrays.
void FinishPolygon()
Decompose polygon in a convex outscribed part and a list of daughter polygons that have to be subtrac...
void SetConvex(Bool_t flag=kTRUE)
Definition TGeoPolygon.h:66
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),...
static Bool_t IsSameWithinTolerance(Double_t a, Double_t b)
Check if two numbers differ with less than a tolerance.
An array of TObjects.
Definition TObjArray.h:37
Int_t GetEntriesFast() const
Definition TObjArray.h:64
Definition TObjArray.h:74
TObject * UncheckedAt(Int_t i) const
Definition TObjArray.h:90
virtual void Delete(Option_t *option="")
Remove all objects from the array AND delete all heap based objects.
TObject * At(Int_t idx) const
Definition TObjArray.h:166
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition TObject.cxx:696
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:893
virtual void Fatal(const char *method, const char *msgfmt,...) const
Issue fatal error message.
Definition TObject.cxx:921
virtual void DrawPolygon(const TGeoPolygon *poly)=0
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17
Double_t Sqrt(Double_t x)
Definition TMath.h:691
Short_t Abs(Short_t d)
Definition TMathBase.h:120