geometry2d.py
1     
2    from typing import TypeVar,List,Sequence; 
3    from math import pi,sqrt,atan2; 
4    from enum import Enum; 
5    #------------------------------------------------------------------------------# 
6     
7    Real=TypeVar('Real',int,float); 
8    #------------------------------------------------------------------------------# 
9     
10   twicePi=2.0*pi; 
11   halfPi=0.5*pi; 
12    
13   r2d=180.0/pi; 
14   d2r=pi/180.0; 
15   #------------------------------------------------------------------------------# 
16    
17   # The constants have to be adapted to the problem 
18   maxLength=10000.0; # area extension 
19   epsilon=0.005 
20   deltaAngles=epsilon/maxLength; # radian 
21   deltaCoordinates=epsilon; 
22    
23   def areCoordinatesEqual(c1,c2:Real)->bool: 
24       return abs(c2-c1)<deltaCoordinates; 
25    
26   def areLengthsEqual(l1,l2:Real)->bool: 
27       return abs(l2-l1)<deltaCoordinates; 
28    
29   def areAnglesEqual(a1,a2:Real)->bool: 
30       return abs(a2-a1)<deltaAngles; 
31   #------------------------------------------------------------------------------# 
32    
33   class BadError(Exception): 
34       def __init__(self,data): 
35           self.Data=data; 
36    
37   class InvalidDataError(Exception): 
38       def __init__(self,data): 
39           self.Data=data; 
40    
41   class InfinityPoints(Exception): 
42       def __init__(self): pass; 
43        
44   class NoPoint(Exception): 
45       def __init__(self): pass; 
46   #------------------------------------------------------------------------------# 
47    
48   class IntersectionKind(Enum): 
49       onePoint=1; 
50       noPoint=2; 
51       infinityPoints=3; 
52   # end IntersectionKind 
53   #------------------------------------------------------------------------------# 
54    
55   class Vectors(object): 
56       def __init__(self,x:Real=0.0,y:Real=0.0)->None: 
57           self.X=x; self.Y=y; 
58    
59       def toClone(self)->'Vectors': 
60           return Vectors(self.X,self.Y); 
61    
62       def __str__(self)->str: 
63           return str(self.X)+', '+str(self.Y); 
64    
65       def __eq__(self,vc:'Vectors')->bool: 
66           """ Test for equality og two vectors  
67           """ 
68           return (areCoordinatesEqual(self.X,vc.X) and 
69                   areCoordinatesEqual(self.Y,vc.Y)); 
70    
71       def lengthOf(self)->Real: 
72           """ Returns the length of the vector 
73           """ 
74           return sqrt(self.X*self.X+self.Y*self.Y); 
75    
76       def toUnitVector(self)->'Vectors': 
77           """ Returns a unit vector in the direction of self  
78           """ 
79           lg=sqrt(self.X*self.X+self.Y*self.Y); 
80           return Vectors(self.X/lg,self.Y/lg); 
81       # end toUnitVector 
82    
83       def __add__(self,vc:'Vectors')->'Vectors': 
84           """ Addition of two vectors  
85           """ 
86           return Vectors(self.X+vc.X,self.Y+vc.Y); 
87    
88       def __sub__(self,vc:'Vectors')->'Vectors': 
89           """ Subtraction of two vectors  
90           """ 
91           return Vectors(self.X-vc.X,self.Y-vc.Y); 
92    
93       def __mul__(self,sc:Real)->'Vectors': 
94           """ Multiplication with scalar to vector  
95           """ 
96           return Vectors(self.X*sc,self.Y*sc);  
97    
98       def smul(self,sc:Real)->None: 
99           """ Multiplication with scalar to self  
100          """ 
101          self.X*=sc; self.Y*=sc;  
102   
103      def __rmul__(self,sc:Real)->'Vectors': 
104          """ Right multiplication with scalar to vector  
105          """ 
106          return Vectors(self.X*sc,self.Y*sc); 
107   
108      def __truediv__(self,sc:Real)->'Vectors': 
109          """ Division with scalar to vector  
110          """ 
111          return Vectors(self.X/sc,self.Y/sc); 
112   
113      def sdiv(self,sc:Real)->None: 
114          """ Division with scalar to self  
115          """ 
116          self.X/=sc; self.Y/=sc; 
117   
118      def dotOp(self,vc:'Vectors')->Real: 
119          """ Dot product of two vectors 
120          """ 
121          return self.X*vc.X+self.Y*vc.Y; 
122   
123      def crossOp(self,vc:'Vectors')->Real: 
124          """ Cross product of two vectors 
125          """ 
126          return self.X*vc.Y-self.Y*vc.X; 
127   
128      def dirAngleOf(self)->Real: 
129          """ Direction angle of the vector relative to the x axis 
130              0<=dirAngleOf<2*Pi 
131          """ 
132          angle=atan2(self.Y,self.X); 
133          if angle<0.0: 
134              return angle+twicePi; 
135          else: 
136              return angle; 
137      # end directionAngle 
138   
139      def sameDirectionOf(self,withV:'Vectors')->bool: 
140          """ Test for equality of the directions of two vectors 
141          """ 
142          angleDiff=abs(withV.dirAngleOf()-self.dirAngleOf()); 
143          if angleDiff<pi: 
144              return angleDiff<deltaAngles; 
145          else: 
146              return (twicePi-angleDiff)<deltaAngles; 
147      # end sameDirectionOf 
148   
149      def isOrthogonalTo(self,vc:'Vectors')->bool: 
150          return abs(self.dotOp(vc))<epsilon; 
151      # end isOrthogonalTo 
152  # end Vectors 
153  #------------------------------------------------------------------------------# 
154   
155  class Points(object): 
156      def __init__(self,x:Real=0.0,y:Real=0.0)->None: 
157          self.X=x; self.Y=y; 
158   
159      def toClone(self)->'Points': 
160          return Points(self.X,self.Y); 
161   
162      def updateOf(self,withP:'Points')->None: 
163          self.X=withP.X; 
164          self.Y=withP.Y; 
165      # end updateOf 
166   
167      def __str__(self)->str: 
168          return str(self.X)+', '+str(self.Y); 
169   
170      def __eq__(self,pt:'Points')->bool: 
171          """ Test for equality of two points 
172          """ 
173          return (areCoordinatesEqual(self.X,pt.X) and 
174                  areCoordinatesEqual(self.Y,pt.Y)); 
175   
176      def __add__(self,vc:Vectors)->'Points': 
177          """ Addition: Point with vector to point 
178          """ 
179          return Points(self.X+vc.X,self.Y+vc.Y); 
180   
181      def __sub__(self,pt:'Points')->'Vectors': 
182          """ Subtraction: Point with point to vector 
183          """ 
184          return Vectors(self.X-pt.X,self.Y-pt.Y); 
185   
186      def distanceOf(self,toPoint:'Points')->Real: 
187          """ Distance between two points 
188          """ 
189          return Vectors.lengthOf(toPoint-self);  
190   
191      def areaOf(self,p1:'Points',p2:'Points')->Real: 
192          """ Area of the triangle with three points 
193          """ 
194          return 0.5*Vectors.crossOp(p1-self,p2-self); 
195  #end Points 
196  #------------------------------------------------------------------------------# 
197   
198  class Lines(object): 
199      def __init__(self,a:Real,b:Real,c:Real)->None: 
200          """ Implicite line equation is A*x+B*y=C 
201              with sqrt(A*A+B*B)=1 
202          """ 
203          d=sqrt(a**2+b**2); 
204          self.A=a/d; self.B=b/d; self.C=c/d; 
205          if areCoordinatesEqual(self.A,0.0): 
206              xlp=0.0; 
207              ylp=self.C/self.B; 
208              xdir=1.0; 
209              ydir=0.0; 
210          elif areCoordinatesEqual(self.B,0.0): 
211              xlp=self.C/self.A; 
212              ylp=0.0; 
213              xdir=0.0; 
214              ydir=1.0; 
215          else: 
216              xlp=0.0; 
217              ylp=self.C/self.B; 
218              xdir= self.B; 
219              ydir= -self.A; 
220          # end if 
221          self.LP=Points(xlp,ylp); 
222          self.DV=Vectors(xdir,ydir); # Unit vector 
223          # end __init__ 
224   
225      def toClone(self)->'Lines': 
226          return Lines(self.A,self.B,self.C); 
227   
228      def __str__(self)->str: 
229          return ('ABC: '+str(self.A)+' '+str(self.B)+' '+str(self.C)+ 
230                  '\n\tLP: '+str(self.LP)+' | DV: '+str(self.DV)); 
231   
232      def __eq__(self,ln:'Lines')->bool: 
233          # Equality of two lines 
234          return (self.LP==ln.LP) and ((self.DV==ln.DV) or (self.DV==(-1.0*ln.DV))); 
235   
236      @staticmethod 
237      def toLine(p1:Points,p2:Points)->'Lines': 
238          """ A line through two points 
239              :raises InvalidDataError 
240                  if p1==p2 
241          """ 
242          if p1==p2: 
243              raise InvalidDataError("toLine with identical points") 
244          a=p1.Y-p2.Y; 
245          b=p2.X-p1.X; 
246          c=a*p1.X + b*p1.Y; 
247          return Lines(a,b,c); 
248      # end toLine 
249   
250      def distanceOf(self,pt:Points)->Real: 
251          """ Signed distance of a point to the line 
252          """ 
253          return self.A*pt.X + self.B*pt.Y - self.C; 
254   
255      def nearestPointOf(self,fromP:Points)->Points: 
256          """ Nearest point on the line from a point 
257          """ 
258          lP=self.LP; 
259          dV=self.DV; 
260          t=Vectors.dotOp(dV,fromP-lP); 
261          lP=lP+t*dV; 
262          return lP; 
263      # end nearestPointOf 
264   
265      def isPointInside(self,pt:Points)->bool: 
266          """ Is a point on the line? 
267              The mathematical equation for a point pt on line is 
268                  pt.X*self.A+pt.Y*self.B-self.C=0 
269          """ 
270          np=self.nearestPointOf(pt); 
271          return pt==np; 
272      # end isPointInside 
273   
274      def areOnSameSide(self,pt1:'Points',pt2:Points)->bool: 
275          """ Are both points on the same side relative to the line? 
276              :Pre self.isPointInside(pt1)==False 
277              :Pre self.isPointInside(pt2)==False 
278          """ 
279          ds1=self.distanceOf(pt1); 
280          ds2=self.distanceOf(pt2); 
281          return (ds1>0 and ds2>0) or ((ds1<0 and ds2<0)) 
282      # end areOnSameSide 
283   
284      def areParallel(self,withLine:'Lines')->bool: 
285          """ Test whether two lines are parallel 
286              Note: Equal lines are considered to be NOT parallel 
287          """ 
288          if (self.LP==withLine.LP): 
289              return False; 
290          d1=self.DV; d2=withLine.DV; 
291          return (d1==d2) or (d1==(-1)*d2); 
292      # end areParallel 
293   
294      def intersectionOf(self,withLine:'Lines')->Points: 
295          """ Returns the point of intersection of two lines 
296              :raises InfinityPoints 
297                  if self==withLine 
298              :raises NoPoint 
299                  if self.areParallel(withLine)==true 
300          """ 
301          L1=self; L2=withLine; 
302          if L1==L2: 
303              raise InfinityPoints; 
304          if L1.areParallel(L2): 
305              raise NoPoint; 
306          lP1=L1.LP.toClone(); ### ??? 
307          lP2=L2.LP; 
308          dV1=L1.DV; 
309          dV2Ortho=Vectors(-L2.DV.Y,L2.DV.X); 
310          q=Vectors.dotOp(dV1,dV2Ortho); 
311          ts=Vectors.dotOp(dV2Ortho,lP2-lP1); 
312          ts/=q; 
313          return lP1+ts*dV1; 
314      # end intersectionOf 
315   
316      def getNormalVector(self)->Vectors: 
317          """ Returns a unit vector perpendicular to the line 
318          """ 
319          return Vectors(self.A,self.B); 
320      # end getNormalVector 
321   
322      def getParallelLine(self,withDistance:Real)->'Lines': 
323          """ Returns the line parallel to self with withDistance 
324          """ 
325          return Lines(self.A,self.B,self.C+withDistance); 
326      # end getParallelLine 
327  # end Lines 
328  #------------------------------------------------------------------------------# 
329   
330  class Rays(object): 
331      def __init__(self,startPoint:Points, 
332                        direction:Vectors)->None: 
333          self.StartPoint=startPoint; 
334          self.Direction=direction.toUnitVector(); 
335      # end __init__ 
336   
337      def toClone(self)->'Rays': 
338          return Rays(self.StartPoint,self.Direction); 
339   
340      def __str__(self)->str: 
341          return 'SP: '+str(self.StartPoint)+' | DV: '+str(self.Direction); 
342   
343      def __eq__(self,withRay:'Rays')->bool: 
344          # Equality of two rays 
345          return ((self.StartPoint==withRay.StartPoint) and 
346                   self.Direction==withRay.Direction); 
347   
348      @staticmethod 
349      def toRay(pt1:Points,pt2:Points)->'Rays': 
350          """ Returns a ray from point pt1 through point pt2 
351          """ 
352          return Rays(pt1,pt2-pt1); 
353       
354      def toLine(self)->Lines: 
355          """ To extend a ray to a line 
356          """ 
357          return (Lines.toLine(self.StartPoint, 
358                  self.StartPoint+self.Direction)); 
359   
360      def nearestPointOf(self,fromP:Points)->Points: 
361          """ Nearest point on the ray from a point 
362          """ 
363          lP=self.StartPoint.toClone(); 
364          dV=self.Direction; 
365          t=Vectors.dotOp(dV,fromP-lP); 
366          if t>=0.0:       
367              lP=lP+t*dV; 
368          return lP; 
369      # end nearestPointOf 
370   
371      def isPointInside(self,pt:Points)->bool: 
372          """ Is a point inside the ray 
373              (Inside means different from the start point) 
374          """ 
375          np=self.nearestPointOf(pt); 
376          if np==self.StartPoint: 
377              return False; 
378          else: 
379              return np==pt; 
380      # end isPointInside 
381   
382      def intersectionOf(self,withRay:'Rays')->Points: 
383          """ Returns the point of intersection of two rays 
384              :raises NoPoint 
385              :raises InfinityPoints 
386          """ 
387          ry1=self; 
388          ry2=withRay; 
389          ln1=ry1.toLine(); 
390          ln2=ry2.toLine(); 
391          if Lines.areParallel(ln1,ln2): 
392              raise NoPoint; 
393          if ln1==ln2: 
394              if (ry1.isPointInside(ry2.StartPoint) or 
395                  ry2.isPointInside(ry1.StartPoint)): 
396                  raise InfinityPoints; 
397              elif ((ry1.StartPoint==ry2.StartPoint) and 
398                    (ry1.Direction==ry2.Direction)): 
399                  raise InfinityPoints; 
400              else: 
401                  raise NoPoint; 
402          # end if 
403          ip=Lines.intersectionOf(ln1,ln2); 
404          if (ry1.isPointInside(ip) and 
405              ry2.isPointInside(ip)): 
406              return ip; 
407          else: 
408              raise NoPoint; 
409          # end if 
410      # end intersectionOf 
411  # end Rays 
412  #------------------------------------------------------------------------------# 
413   
414  class Segments(object): 
415      def __init__(self,p1:Points,p2:Points)->None: 
416          """ The oriented line segment from point p1 to point p2 
417          """ 
418          self.P1=p1; self.P2=p2; 
419   
420      def toClone(self)->'Segments': 
421          return Segments(self.P1,self.P2); 
422   
423      def update(self,withSegment:'Segments')->None: 
424          self.P1=withSegment.P1; 
425          self.P2=withSegment.P2; 
426      # end update 
427   
428      def __eq__(self,sgm:'Segments')->bool: 
429          """ Test for oriented equality of two segments 
430          """ 
431          return (self.P1==sgm.P1) and (self.P2==sgm.P2); 
432   
433      def __str__(self)->str: 
434          return str(self.P1)+' | '+str(self.P2); 
435   
436      def toLine(self)->Lines: 
437          """ Returns the extension of the segment to a line 
438          """ 
439          return Lines.toLine(self.P1,self.P2); 
440   
441      def isPointInside(self,pt:Points)->bool: 
442          """ Is a point inside the segment? 
443          """ 
444          ps1=self.P1; 
445          ps2=self.P2; 
446          rs1=Rays.toRay(ps1,ps2); 
447          rs2=Rays.toRay(ps2,ps1); 
448          return rs1.isPointInside(pt) and rs2.isPointInside(pt); 
449      # end isPointInside 
450   
451      def distanceOf(self,pt:Points)->bool: 
452          """ Returns the signed distance of a point from the segment 
453          :Pre pt is not a point of the segment 
454          """ 
455          From=pt; To=self; 
456          ln=To.toLine(); 
457          dOrtho=Vectors(ln.A,ln.B); 
458          if (Vectors.dotOp(From-To.P2,To.P2-To.P1)>0.0): 
459              sgn=1 if (Vectors.dotOp(dOrtho,From-To.P2)>=0) else -1; 
460              return sgn*From.distanceOf(To.P2); 
461          elif (Vectors.dotOp(From-To.P1,To.P1-To.P2)>0.0): 
462              sgn=1 if (Vectors.dotOp(dOrtho,From-To.P1)>=0) else -1; 
463              return sgn*From.distanceOf(To.P1); 
464          else: 
465              return ln.distanceOf(From); 
466          # end if 
467      # end distanceOf 
468   
469      def intersectionOf_withRay(self,withRay:Rays)->Points: 
470          """ Returns the point of intersection of the segment with the ray 
471              :raises NoPoint 
472              :raises InfinityPoints 
473          """ 
474          ps1=self.P1; 
475          ps2=self.P2; 
476          rs1=Rays.toRay(ps1,ps2); 
477          rs2=Rays.toRay(ps2,ps1); 
478   
479          intersectionKind1=IntersectionKind.onePoint; 
480          intersectionKind2=IntersectionKind.onePoint; 
481   
482          try: 
483              ip1=withRay.intersectionOf(rs1); 
484              if (ip1==ps1) or (ip1==ps2): 
485                  raise NoPoint; 
486          except NoPoint: 
487              intersectionKind1=IntersectionKind.noPoint; 
488          except InfinityPoints: 
489              intersectionKind1=IntersectionKind.infinityPoints; 
490          # end try 
491   
492          try: 
493              ip2=withRay.intersectionOf(rs2); 
494              if (ip2==ps1) or (ip2==ps2): 
495                  raise NoPoint; 
496          except NoPoint: 
497              intersectionKind2=IntersectionKind.noPoint; 
498          except InfinityPoints: 
499              intersectionKind1=IntersectionKind.infinityPoints; 
500          # end try 
501   
502          if ((intersectionKind1==IntersectionKind.onePoint) and 
503              (intersectionKind2==IntersectionKind.onePoint)): 
504              assert(ip1==ip2); # TBD 
505              return ip1; 
506          elif ((intersectionKind1==IntersectionKind.noPoint) and 
507                (intersectionKind2==IntersectionKind.noPoint)): 
508              raise NoPoint; 
509          elif ((intersectionKind1==IntersectionKind.noPoint) and 
510                (intersectionKind2==IntersectionKind.onePoint)): 
511              raise NoPoint; 
512          elif ((intersectionKind2==IntersectionKind.noPoint) and 
513                (intersectionKind1==IntersectionKind.onePoint)): 
514              raise NoPoint; 
515          else: 
516              raise InfinityPoints; # TBD 
517          # end if 
518      # end intersectionOf_withRay 
519   
520      def intersectionOf(self,withSegment:'Segments')->Points: 
521          """ Returns the point of intersection of two segments 
522              :raises NoPoint 
523              :raises InfinityPoints 
524          """ 
525          ps1=self.P1; 
526          ps2=self.P2; 
527          rs1=Rays.toRay(ps1,ps2); 
528          rs2=Rays.toRay(ps2,ps1); 
529          try: 
530              ip1=withSegment.intersectionOf_withRay(rs1); 
531              ip2=withSegment.intersectionOf_withRay(rs2); 
532          except InfinityPoints: 
533              if (withSegment.isPointInside(ps1) or 
534                  withSegment.isPointInside(ps2)): 
535                  raise InfinityPoints; 
536              else: 
537                  raise NoPoint; 
538              # end if 
539          except NoPoint: 
540              raise NoPoint; 
541          else: 
542              assert ip1==ip2; 
543              return ip1; 
544          # end try 
545      # end intersectionOf 
546   
547      def bisectorOf(self)->Lines: 
548          """ Returns the bisector to self, 
549              i.e. the line that is perpendicular 
550              to self and goes through its middle 
551          """ 
552          ln=self.toLine(); 
553          xMid=(self.P1.X+self.P2.X)/2.0; 
554          yMid=(self.P1.Y+self.P2.Y)/2.0; 
555          return Lines(-ln.B, 
556                        ln.A, 
557                       -ln.B*xMid + ln.A*yMid); 
558      # end bisectorOf 
559   
560      @staticmethod 
561      def areNotIntersecting(segments:Sequence['Segments'])->None: 
562          """ The segments in the list must not intersect 
563              :raises InvalidDataError 
564                  if two segments are intersecting 
565          """ 
566          for sgo in segments: 
567              for sgi in segments: 
568                  if not sgo==sgi: 
569                      try: 
570                          sgo.intersectionOf(sgi); 
571                      except NoPoint: 
572                          pass; # Requirement 
573                      except InfinityPoints: 
574                          raise InvalidDataError('Overlaping segments'); 
575                      else: 
576                          raise InvalidDataError('Overlaping segments'); 
577                      # end try 
578                  # end if 
579              # end for 
580          # end for 
581      # end areNotIntersecting 
582  # end Segments 
583  #------------------------------------------------------------------------------# 
584   
585  class Sectors(object): 
586      def __init__(self,startPoint:Points, 
587                        direction1:Vectors, 
588                        direction2:Vectors)->None: 
589          self.StartPoint=startPoint; 
590          self.Direction1=direction1.toUnitVector(); 
591          self.Direction2=direction2.toUnitVector(); 
592      # end __init__ 
593   
594      def toClone(self)->'Sectors': 
595          return Sectors(self.StartPoint,self.Direction1,self.Direction2); 
596   
597      def __str__(self)->str: 
598          return ('SP: '+str(self.StartPoint)+ 
599                  ' | DV1: '+str(self.Direction1)+ 
600                  ' | DV2: '+str(self.Direction2)); 
601   
602      def __eq__(self,withSector:'Sectors')->bool: 
603          return ((self.StartPoint==withSector.StartPoint) and 
604                   ((self.Direction1==withSector.Direction1) and 
605                   (self.Direction2==withSector.Direction2)) or 
606                   ((self.Direction1==withSector.Direction2) and 
607                   (self.Direction2==withSector.Direction1))); 
608   
609      @staticmethod 
610      def toSector(pt0:Points,pt1:Points,pt2:Points)->'Sectors': 
611          return Sectors(pt0,pt1-pt0,pt2-pt0); 
612   
613      def openingAngleOf(self)->Real: 
614          """ Gibt den nicht-orientierten Öffnungswinkel zurück 
615          """ 
616          return abs(self.Direction2.dirAngleOf()-self.Direction1.dirAngleOf()); 
617  # end Sectors 
618  #------------------------------------------------------------------------------# 
619   
620  class SectorWithSegments(Sectors): 
621      def __init__(self,startPoint:Points, 
622                        direction1:Vectors, 
623                        direction2:Vectors)->None: 
624          self.Segments=[]; 
625          Sectors.__init__(self,startPoint,direction1,direction2); 
626      # end __init__ 
627   
628      def __str__(self)->str: 
629          output=''; 
630          output+=Sectors.__str__(self); 
631          output+=' || '; 
632          for sgm in self.Segments: 
633              output+=str(sgm)+' | '; 
634          return(output); 
635      # end __str__ 
636   
637      def add(self,segment:Segments)->bool: 
638          """  Tries to add the segment to self 
639          """ 
640          p0=self.StartPoint; 
641          p1=segment.P1; p2=segment.P2; 
642          sec=Sectors.toSector(p0,p1,p2); 
643          if self==sec: 
644              self.Segments.append(segment); 
645              return True; 
646          else: 
647              return False; 
648          # end if 
649      # end add 
650   
651      @staticmethod 
652      def addToList(secwsgm:Sequence['SectorWithSegments'], 
653                    sgm:Segments)->None: 
654          """ Adds the segment to one of the sectors of the list 
655              :raises InvalidDataError 
656                  if there is no sector suitable for the segment 
657          """ 
658          added=False; 
659          for sec in secwsgm: 
660              added=False; 
661              if sec.add(sgm): 
662                  added=True; 
663                  break; 
664          # end for 
665          if not added: 
666              raise InvalidDataError('addToList'); 
667      # end addToList 
668   
669      def getMinimalDistanceSegment(self)->Segments: 
670          """ Returns the segment that is nearest to the StartPoint 
671              of the sector 
672          """ 
673          segments=self.Segments; 
674          toPoint=self.StartPoint; 
675          minSegment=segments[0]; 
676          minlg1=Vectors.lengthOf(minSegment.P1-toPoint); 
677          minlg2=Vectors.lengthOf(minSegment.P2-toPoint); 
678   
679          for sgm in segments: 
680              lg1=Vectors.lengthOf(sgm.P1-toPoint); 
681              lg2=Vectors.lengthOf(sgm.P2-toPoint); 
682   
683              if ((lg1<minlg1 and lg2<=minlg2) or 
684                  (lg1<=minlg1 and lg2<minlg2)): 
685                  minlg1=lg1; 
686                  minlg2=lg2; 
687                  minSegment.update(sgm); 
688              # end if 
689          # end for 
690   
691          return minSegment; 
692      # end getMinimalDistanceSegment 
693  # end SectorWithSegments 
694  #------------------------------------------------------------------------------# 
695   
696  class Triangles(object): 
697      def __init__(self,p0:Points,p1:Points,p2:Points)->None: 
698          self.P0=p0; self.P1=p1; self.P2=p2; 
699   
700      def __str__(self)->str: 
701          return str(self.P0)+' | '+str(self.P1)+' | '+str(self.P2); 
702   
703      def areaOf(self)->Real: 
704          return Points.areaOf(self.P0,self.P1,self.P2); 
705   
706      def isPointInside(self,pt:Points)->bool: 
707          areaT=Points.areaOf(self.P0,self.P1,self.P2); 
708          # Barycentric coordinates with (alpha+beta+gamma)=1 
709          alpha=Points.areaOf(pt,self.P1,self.P2)/areaT; 
710          beta =Points.areaOf(pt,self.P2,self.P0)/areaT; 
711          gamma=Points.areaOf(pt,self.P0,self.P1)/areaT; 
712           
713          return (0.0<alpha<1.0) and (0.0<beta<1.0) and (0.0<gamma<1.0); 
714              #! No use of some delta - TBD 
715      # end isPointInside 
716  # end Triangles 
717  #------------------------------------------------------------------------------# 
718   
719  class Polygons(object): 
720      def __init__(self,points:Sequence[Points])->None: 
721          self.NbPoints=len(points); 
722          self.Points=[]; 
723          self.Points+=points; 
724      # end __init__ 
725   
726      def toClone(self)->'Polygons': 
727          return Polygons(self.Points); 
728   
729      def __str__(self)->str: 
730          output=''; 
731          for pt in self.Points: 
732              output+=str(pt)+' | '; 
733   
734          return(output);           
735      # end __str__ 
736   
737      def distanceOf(self,pt:Points)->Real: 
738          dist=-1; 
739          for Ix in range(self.NbPoints-1): 
740              p1=self.Points[Ix]; 
741              p2=self.Points[Ix+1]; 
742              sgm=Segments(p1,p2); 
743              dist=min(maxLength,sgm.distanceOf(pt));            
744          # end for 
745          p1=self.Points[1]; 
746          p2=self.Points[self.NbPoints-1]; 
747          sgm=Segments(p1,p2); 
748          dist=min(dist,sgm.distanceOf(pt)); 
749          return dist; 
750      # end distanceOf 
751   
752      def areaOf(self)->Real: 
753          """ Returns the signed area of the polygon 
754          """ 
755          p0=self.Points[0]; 
756          area=0.0; 
757          for Ix in range(1,self.NbPoints-1): 
758              p1=self.Points[Ix]; 
759              p2=self.Points[Ix+1]; 
760              area+=Vectors.crossOp(p1-p0,p2-p0); 
761          # end for 
762          return 0.5*area; 
763      # end areaOf 
764   
765      def isPointOnside (self,pt0:Points)->bool: 
766          """ Is a point onside the polygon? 
767              Onside means: The point is inside or on the edges of the polygon 
768              This algorithm is adapted from 
769              'The GNAT Components Collection', www.adacore.com 
770          """ 
771          Ie=self.NbPoints-1; 
772          onInside= False; 
773          deltaY=-1.0; 
774          points=self.Points; 
775   
776          for Ix in range(Ie): 
777              deltaY=pt0.Y-points[Ix].Y; 
778   
779              if (((0.0<=deltaY and pt0.Y<points[Ie].Y) 
780                    or (points[Ie].Y<=pt0.Y and deltaY<0.0)) 
781                  and (pt0.X-points[Ix].X<(points[Ie].X-points[Ix].X)*deltaY 
782                    /(points[Ie].Y-points[Ix].Y))): 
783   
784                  onInside=not onInside; 
785           # end if 
786        # end for 
787   
788          return onInside; 
789     # end isPointOnside; 
790  # end Polygons 
791  #------------------------------------------------------------------------------# 
792   
793  class Circles(object): 
794      class NoCircle(Exception): 
795          def __init__(self): pass;  
796   
797      def __init__(self,point:Points,radius:Real)->None: 
798          self.C=point; 
799          self.R=radius; 
800      # end def 
801   
802      def toClone(self)->'Circles': 
803          return Circles(self.C,self.R); 
804   
805      def __str__(self)->str: 
806          return str(self.C)+' | '+str(self.R); 
807   
808      @staticmethod 
809      def toCircle(pt1:Points,pt2:Points,pt3:Points)->'Circles': 
810          """ Creates the circle through 3 points 
811              :raises NoCircle 
812                  if points are on a line 
813          """ 
814          sgm1=Segments(pt1,pt2); 
815          sgm2=Segments(pt2,pt3); 
816          bis1=sgm1.bisectorOf(); 
817          bis2=sgm2.bisectorOf(); 
818          try: 
819              center=Lines.intersectionOf(bis1,bis2); 
820          except (NoPoint,InfinityPoints): 
821              raise Circles.NoCircle; 
822          # end try 
823          radius=pt1.distanceOf(center); 
824          return Circles(center,radius); 
825      # end toCircle 
826   
827      def areaOf(self)->Real: 
828          return pi*self.R*self.R; 
829   
830      def distanceOf(self,pt:Points)->Real: 
831          """ Returns the signed distance of a point to the circle 
832          """ 
833          d=Points.distanceOf(self.C,pt); 
834          return d-self.R; 
835      # end distanceOf 
836   
837      def isPointInside(self,pt:Points)->bool: 
838          dist=self.distanceOf(pt); 
839          return (dist<0) and not areLengthsEqual(dist,0); 
840      # end isPointInside 
841  # end Circles 
842  #------------------------------------------------------------------------------# 
843   
844  class PolygonTrains(object): # Polygonzug 
845   
846      def __init__(self,points:Sequence[Points])->None: 
847          self.Points=[]; 
848          self.Points+=points; 
849      # end __init__ 
850   
851      def getClone(self)->'PolygonTrains': 
852          return PolygonTrains(self.Points); 
853   
854      def __str__(self)->str: 
855          output=''; 
856          for pt in self.Points: 
857              output+=str(pt)+' | '; 
858          return(output); 
859      # end __str__ 
860   
861      def getSegments(self)->List[Segments]: 
862          """ Returns the segments of the polygon train 
863          """ 
864          la=self.Points; 
865          fIx=0; lIx=len(la)-1; 
866          lf=la[fIx:lIx]; ll=la[fIx+1:lIx+1]; 
867          return [Segments(*tp) for tp in list(zip(lf,ll))]; 
868      # end getSegments 
869  # end PolygonTrains 
870  #------------------------------------------------------------------------------# 
871   
872  class ConesOfRays(object): # Strahlenbueschel 
873   
874      def __init__(self,startPoint:Points)->None: 
875          self.StartPoint=startPoint; 
876          self.Rays=[]; 
877      # end __init__ 
878   
879      def __str__(self)->str: 
880          output=''; 
881          for ray in self.Rays: 
882              angle=ray.Direction.dirAngleOf()*r2d; 
883              output+='\n'+str(angle)+' - '+str(ray)+' | '; 
884          # end for 
885          return(output); 
886      # end __str__ 
887   
888      def add(self,ray:Rays)->None: 
889          """ Add the ray to self 
890              :raises InvalidDataError 
891                  if start points of rays are different 
892          """ 
893          if not self.StartPoint==ray.StartPoint: 
894              raise InvalidDataError('different start points'); 
895          for elem in self.Rays: 
896              if elem==ray: return None; # No dublicates 
897          self.Rays.insert(0,ray); 
898      # end add 
899   
900      @staticmethod 
901      def createFromPoints(startPoint:Points, 
902                           points:Sequence[Points])->'ConesOfRays': 
903          """ Creates a cone of rays using a start point 
904              for the rays and some other points 
905          """ 
906          cone=ConesOfRays(startPoint); 
907          for pt in points: 
908              dirc=(pt-startPoint).toUnitVector(); 
909              cone.add(Rays(startPoint,dirc)); 
910          # end for 
911   
912          return cone; 
913      # end createFromPoints 
914  # ConesOfRays 
915  #------------------------------------------------------------------------------# 
916   
917      def divideIntersectedSegments(self,sgms:List[Segments])->None: 
918          """ Schneidet ein Strahl ein Segment im Innern, 
919              so wird dieses Segment in zwei Einzelsegmente zerlegt 
920          """ 
921          for ray in self.Rays: 
922              for sgm in sgms[:]: 
923                  try: 
924                      ip=sgm.intersectionOf_withRay(ray); 
925                  except NoPoint: 
926                      pass; 
927                  except InfinityPoints: 
928                      sgms.remove(sgm); 
929                      # Segment is part of ray; no extension, no visibility 
930                  else: 
931                      sgm1=Segments(sgm.P1,ip); 
932                      sgms.append(sgm1); 
933                      sgm2=Segments(sgm.P2,ip); 
934                      sgms.append(sgm2); 
935                      sgms.remove(sgm); 
936                  # end try 
937              # end for 
938          # end for 
939      # divideIntersectedSegments 
940  # end ConesOfRays 
941  #------------------------------------------------------------------------------# 
942  #------------------------------------------------------------------------------# 
943   
944  if __name__ == '__main__': 
945      print('\n# Beginning main ...'); 
946   
947      print('\n# Finished main.\n'); 
948  # end if main 
949   
950  #------------------------------------------------------------------------------# 
951  #------------------------------------------------------------------------------#