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 Shapes_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)
105 delete[] fInd;
106 if (fIndc)
107 delete[] fIndc;
108 if (fDaughters) {
110 delete fDaughters;
111 }
112}
113////////////////////////////////////////////////////////////////////////////////
114/// Computes area of the polygon in [length^2].
115
117{
118 Int_t ic, i, j;
119 Double_t area = 0;
120 // Compute area of the convex part
121 for (ic = 0; ic < fNvert; ic++) {
122 i = fInd[ic];
123 j = fInd[(ic + 1) % fNvert];
124 area += 0.5 * (fX[i] * fY[j] - fX[j] * fY[i]);
125 }
126 return TMath::Abs(area);
127}
128
129////////////////////////////////////////////////////////////////////////////////
130/// Check if a point given by X = point[0], Y = point[1] is inside the polygon.
131
133{
134 Int_t i;
135 TGeoPolygon *poly;
136 for (i = 0; i < fNconvex; i++)
137 if (!IsRightSided(point, fIndc[i], fIndc[(i + 1) % fNconvex]))
138 return kFALSE;
139 if (!fDaughters)
140 return kTRUE;
142 for (i = 0; i < nd; i++) {
143 poly = (TGeoPolygon *)fDaughters->UncheckedAt(i);
144 if (poly->Contains(point))
145 return kFALSE;
146 }
147 return kTRUE;
148}
149
150////////////////////////////////////////////////////////////////////////////////
151/// Check polygon convexity.
152
154{
155 if (fNvert == 3) {
156 SetConvex();
157 return;
158 }
159 Int_t j, k;
160 Double_t point[2];
161 for (Int_t i = 0; i < fNvert; i++) {
162 j = (i + 1) % fNvert;
163 k = (i + 2) % fNvert;
164 point[0] = fX[fInd[k]];
165 point[1] = fY[fInd[k]];
166 if (!IsRightSided(point, fInd[i], fInd[j]))
167 return;
168 }
169 SetConvex();
170}
171
172////////////////////////////////////////////////////////////////////////////////
173/// Draw the polygon.
174
176{
177 if (!gGeoManager)
178 return;
180}
181
182////////////////////////////////////////////////////////////////////////////////
183/// Decompose polygon in a convex outscribed part and a list of daughter
184/// polygons that have to be subtracted to get the actual one.
185
187{
189 // check convexity
190 ConvexCheck();
191 // find outscribed convex polygon indices
193 if (IsConvex()) {
194 // printf(" -> polygon convex -> same indices\n");
195 memcpy(fIndc, fInd, fNvert * sizeof(Int_t));
196 return;
197 }
198 // printf(" -> polygon NOT convex\n");
199 // make daughters if necessary
200 if (IsConvex())
201 return;
202 // ... algorithm here
203 if (!fDaughters)
204 fDaughters = new TObjArray();
205 TGeoPolygon *poly = nullptr;
206 Int_t indconv = 0;
207 Int_t indnext, indback;
208 Int_t nskip;
209 while (indconv < fNconvex) {
210 indnext = (indconv + 1) % fNconvex;
211 nskip = fIndc[indnext] - fIndc[indconv];
212 if (nskip < 0)
213 nskip += fNvert;
214 if (nskip == 1) {
215 indconv++;
216 continue;
217 }
218 // gap -> make polygon
219 poly = new TGeoPolygon(nskip + 1);
220 poly->SetXY(fX, fY);
221 poly->SetNextIndex(fInd[fIndc[indconv]]);
222 poly->SetNextIndex(fInd[fIndc[indnext]]);
223 indback = fIndc[indnext] - 1;
224 if (indback < 0)
225 indback += fNvert;
226 while (indback != fIndc[indconv]) {
227 poly->SetNextIndex(fInd[indback]);
228 indback--;
229 if (indback < 0)
230 indback += fNvert;
231 }
232 poly->FinishPolygon();
234 indconv++;
235 }
236 for (indconv = 0; indconv < fNconvex; indconv++)
237 fIndc[indconv] = fInd[fIndc[indconv]];
238}
239
240////////////////////////////////////////////////////////////////////////////////
241/// Fill list of vertices into provided arrays.
242
244{
245 memcpy(x, fX, fNvert * sizeof(Double_t));
246 memcpy(y, fY, fNvert * sizeof(Double_t));
247}
248
249////////////////////////////////////////////////////////////////////////////////
250/// Fill list of vertices of the convex outscribed polygon into provided arrays.
251
253{
254 for (Int_t ic = 0; ic < fNconvex; ic++) {
255 x[ic] = fX[fIndc[ic]];
256 y[ic] = fY[fIndc[ic]];
257 }
258}
259
260////////////////////////////////////////////////////////////////////////////////
261/// Check if POINT is right-sided with respect to the segment defined by IND1 and IND2.
262
263Bool_t TGeoPolygon::IsRightSided(const Double_t *point, Int_t ind1, Int_t ind2) const
264{
265 Double_t dot = (point[0] - fX[ind1]) * (fY[ind2] - fY[ind1]) - (point[1] - fY[ind1]) * (fX[ind2] - fX[ind1]);
266 if (!IsClockwise())
267 dot = -dot;
268 if (dot < -1.E-10)
269 return kFALSE;
270 return kTRUE;
271}
272
273////////////////////////////////////////////////////////////////////////////////
274/// Check if a segment [0..fNvert-1] belongs to the outscribed convex pgon.
275
277{
278 if (i2 < 0)
279 i2 = (i1 + 1) % fNvert;
280 Double_t point[2];
281 for (Int_t i = 0; i < fNvert; i++) {
282 if (i == i1 || i == i2)
283 continue;
284 point[0] = fX[fInd[i]];
285 point[1] = fY[fInd[i]];
286 if (!IsRightSided(point, fInd[i1], fInd[i2]))
287 return kFALSE;
288 }
289 return kTRUE;
290}
291
292////////////////////////////////////////////////////////////////////////////////
293/// Check for illegal crossings between non-consecutive segments
294
296{
297 if (fNvert < 4)
298 return kFALSE;
299 Bool_t is_illegal = kFALSE;
300 Double_t x1, y1, x2, y2, x3, y3, x4, y4;
301 for (Int_t i = 0; i < fNvert - 2; i++) {
302 // Check segment i
303 for (Int_t j = i + 2; j < fNvert; j++) {
304 // Versus segment j
305 if (i == 0 && j == (fNvert - 1))
306 continue;
307 x1 = fX[i];
308 y1 = fY[i];
309 x2 = fX[i + 1];
310 y2 = fY[i + 1];
311 x3 = fX[j];
312 y3 = fY[j];
313 x4 = fX[(j + 1) % fNvert];
314 y4 = fY[(j + 1) % fNvert];
315 if (TGeoShape::IsSegCrossing(x1, y1, x2, y2, x3, y3, x4, y4)) {
316 Error("IsIllegalCheck", "Illegal crossing of segment %d vs. segment %d", i, j);
317 is_illegal = kTRUE;
318 }
319 }
320 }
321 return is_illegal;
322}
323
324////////////////////////////////////////////////////////////////////////////////
325/// Compute indices for the outscribed convex polygon.
326
328{
329 fNconvex = 0;
330 Int_t iseg = 0;
331 Int_t ivnew;
332 Bool_t conv;
333 Int_t *indconv = new Int_t[fNvert];
334 memset(indconv, 0, fNvert * sizeof(Int_t));
335 while (iseg < fNvert) {
336 if (!IsSegConvex(iseg)) {
337 if (iseg + 2 > fNvert)
338 break;
339 ivnew = (iseg + 2) % fNvert;
340 conv = kFALSE;
341 // check iseg with next vertices
342 while (ivnew) {
343 if (IsSegConvex(iseg, ivnew)) {
344 conv = kTRUE;
345 break;
346 }
347 ivnew = (ivnew + 1) % fNvert;
348 }
349 if (!conv) {
350 // Error("OutscribedConvex","NO convex line connection to vertex %d\n", iseg);
351 iseg++;
352 continue;
353 }
354 } else {
355 ivnew = (iseg + 1) % fNvert;
356 }
357 // segment belonging to convex outscribed polygon
358 if (!fNconvex)
359 indconv[fNconvex++] = iseg;
360 else if (indconv[fNconvex - 1] != iseg)
361 indconv[fNconvex++] = iseg;
362 if (iseg < fNvert - 1)
363 indconv[fNconvex++] = ivnew;
364 if (ivnew < iseg)
365 break;
366 iseg = ivnew;
367 }
368 if (!fNconvex) {
369 delete[] indconv;
370 Fatal("OutscribedConvex", "cannot build outscribed convex");
371 return;
372 }
373 fIndc = new Int_t[fNvert];
374 memcpy(fIndc, indconv, fNconvex * sizeof(Int_t)); // does not contain real indices yet
375 delete[] indconv;
376}
377
378////////////////////////////////////////////////////////////////////////////////
379/// Compute minimum distance from POINT to any segment. Returns segment index.
380
381Double_t TGeoPolygon::Safety(const Double_t *point, Int_t &isegment) const
382{
383 Int_t i1, i2;
384 Double_t p1[2], p2[3];
385 Double_t lsq, ssq, dx, dy, dpx, dpy, u;
386 Double_t safe = 1E30;
387 Int_t isegmin = 0;
388 for (i1 = 0; i1 < fNvert; i1++) {
390 isegment = isegmin;
391 return 0.;
392 }
393 i2 = (i1 + 1) % fNvert;
394 p1[0] = fX[i1];
395 p1[1] = fY[i1];
396 p2[0] = fX[i2];
397 p2[1] = fY[i2];
398
399 dx = p2[0] - p1[0];
400 dy = p2[1] - p1[1];
401 dpx = point[0] - p1[0];
402 dpy = point[1] - p1[1];
403
404 lsq = dx * dx + dy * dy;
406 ssq = dpx * dpx + dpy * dpy;
407 if (ssq < safe) {
408 safe = ssq;
409 isegmin = i1;
410 }
411 continue;
412 }
413 u = (dpx * dx + dpy * dy) / lsq;
414 if (u > 1) {
415 dpx = point[0] - p2[0];
416 dpy = point[1] - p2[1];
417 } else {
418 if (u >= 0) {
419 dpx -= u * dx;
420 dpy -= u * dy;
421 }
422 }
423 ssq = dpx * dpx + dpy * dpy;
424 if (ssq < safe) {
425 safe = ssq;
426 isegmin = i1;
427 }
428 }
429 isegment = isegmin;
430 safe = TMath::Sqrt(safe);
431 // printf("== segment %d: (%f, %f) - (%f, %f) safe=%f\n", isegment,
432 // fX[isegment],fY[isegment],fX[(isegment+1)%fNvert],fY[(isegment+1)%fNvert],safe);
433 return safe;
434}
435
436////////////////////////////////////////////////////////////////////////////////
437/// Sets the next polygone index. If index<0 sets all indices consecutive
438/// in increasing order.
439
441{
442 Int_t i;
443 if (index < 0) {
444 for (i = 0; i < fNvert; i++)
445 fInd[i] = i;
446 return;
447 }
448 if (fNconvex >= fNvert) {
449 Error("SetNextIndex", "all indices already set");
450 return;
451 }
452 fInd[fNconvex++] = index;
453 if (fNconvex == fNvert) {
454 if (!fX || !fY)
455 return;
456 Double_t area = 0.0;
457 for (i = 0; i < fNvert; i++)
458 area += fX[fInd[i]] * fY[fInd[(i + 1) % fNvert]] - fX[fInd[(i + 1) % fNvert]] * fY[fInd[i]];
459 if (area < 0)
461 else
463 }
464}
465
466////////////////////////////////////////////////////////////////////////////////
467/// Set X/Y array pointer for the polygon and daughters.
468
470{
471 Int_t i;
472 fX = x;
473 fY = y;
474 Double_t area = 0.0;
475 for (i = 0; i < fNvert; i++)
476 area += fX[fInd[i]] * fY[fInd[(i + 1) % fNvert]] - fX[fInd[(i + 1) % fNvert]] * fY[fInd[i]];
477 if (area < 0)
479 else
481
482 if (!fDaughters)
483 return;
484 TGeoPolygon *poly;
486 for (i = 0; i < nd; i++) {
487 poly = (TGeoPolygon *)fDaughters->At(i);
488 if (poly)
489 poly->SetXY(x, y);
490 }
491}
constexpr Bool_t kFALSE
Definition RtypesCore.h:101
constexpr Bool_t kTRUE
Definition RtypesCore.h:100
const char Option_t
Definition RtypesCore.h:66
#define ClassImp(name)
Definition Rtypes.h:377
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t index
Option_t Option_t TPoint TPoint const char x2
Option_t Option_t TPoint TPoint const char x1
Option_t Option_t TPoint TPoint const char y2
Option_t Option_t TPoint TPoint const char y1
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:19
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.
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[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:26
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:27
Double_t * fY
pointer to list of current X coordinates of vertices
Definition TGeoPolygon.h:30
void Draw(Option_t *option="") override
Draw the polygon.
Bool_t IsClockwise() const
Definition TGeoPolygon.h:58
void OutscribedConvex()
Compute indices for the outscribed convex polygon.
Double_t * fX
Definition TGeoPolygon.h:29
~TGeoPolygon() override
Destructor.
Int_t * fIndc
Definition TGeoPolygon.h:28
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:31
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),...
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:31
Int_t GetEntriesFast() const
Definition TObjArray.h:58
void Delete(Option_t *option="") override
Remove all objects from the array AND delete all heap based objects.
TObject * At(Int_t idx) const override
Definition TObjArray.h:164
TObject * UncheckedAt(Int_t i) const
Definition TObjArray.h:84
Definition TObjArray.h:68
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition TObject.cxx:780
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:987
virtual void Fatal(const char *method, const char *msgfmt,...) const
Issue fatal error message.
Definition TObject.cxx:1015
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)
Returns the square root of x.
Definition TMath.h:662
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:123