Logo ROOT   6.16/01
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
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 "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;
96}
97
98////////////////////////////////////////////////////////////////////////////////
99/// Destructor
100
102{
103 if (fInd) delete [] fInd;
104 if (fIndc) delete [] fIndc;
105 if (fDaughters) {
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++) {
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
249Bool_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
356Double_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++) {
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}
static double p1(double t, double a, double b)
static double p2(double t, double a, double b, double c)
static const double x2[5]
static const double x4[22]
static const double x1[5]
static const double x3[11]
int Int_t
Definition: RtypesCore.h:41
const Bool_t kFALSE
Definition: RtypesCore.h:88
bool Bool_t
Definition: RtypesCore.h:59
double Double_t
Definition: RtypesCore.h:55
const Bool_t kTRUE
Definition: RtypesCore.h:87
const char Option_t
Definition: RtypesCore.h:62
#define ClassImp(name)
Definition: Rtypes.h:363
R__EXTERN TGeoManager * gGeoManager
Definition: TGeoManager.h:572
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:59
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.
Definition: TGeoPolygon.cxx:64
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].
@ kGeoFinishPolygon
Definition: TGeoPolygon.h:24
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 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:58
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.
Int_t fNvert
Definition: TGeoPolygon.h:29
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:63
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
static Bool_t IsSameWithinTolerance(Double_t a, Double_t b)
Check if two numbers differ with less than a tolerance.
Definition: TGeoShape.cxx:326
An array of TObjects.
Definition: TObjArray.h:37
Int_t GetEntriesFast() const
Definition: TObjArray.h:64
void Add(TObject *obj)
Definition: TObjArray.h:73
TObject * UncheckedAt(Int_t i) const
Definition: TObjArray.h:89
virtual void Delete(Option_t *option="")
Remove all objects from the array AND delete all heap based objects.
Definition: TObjArray.cxx:355
TObject * At(Int_t idx) const
Definition: TObjArray.h:165
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition: TObject.cxx:694
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:880
virtual void Fatal(const char *method, const char *msgfmt,...) const
Issue fatal error message.
Definition: TObject.cxx:908
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:679
Short_t Abs(Short_t d)
Definition: TMathBase.h:120