🚀 go-pugleaf

This is a multi-part message in MIME format.
--------------7EC501A120D126B6D5D2597A
Content-Type: text/plain; charset=us-ascii
Content-Transfer-Encoding: 7bit

Hi,

Well I 've decided to put back into a community that has helped me out
so much over the years.
I've developed a set of really fast geometric routines that I would like
to share with everyone.
Please distribute anyway you like and feel free to use the routines in
your applications
I would also be very happy to hear your comments and ideas or bug
reports (if any :) )


Well I hope I can do this again sometime.


Arash Partow
www.arash.da.ru

--------------7EC501A120D126B6D5D2597A
Content-Type: text/plain; charset=iso-8859-1;
 name="Fastgeo.pas"
Content-Transfer-Encoding: 8bit
Content-Disposition: inline;
 filename="Fastgeo.pas"


(***********************************************************)
(*                                                         *)
(*       Super fast 2D geometric algorithims               *)
(*       Release Ver 1.0  10/10/1999                       *)
(*                                                         *)
(*       Pascal Implementation by Arash Partow             *)
(*                                                         *)
(*       This code is governed by GPL                      *)
(*       You should NOT use this code                      *)
(*       if you are not familiar with                      *)
(*       the term "G.P.L"                                  *)
(*                                                         *)
(*              ("`-''-'").___..--''"`-._                  *)
(*                6_ 6  )   `-.  (    ).`-..,...'`)        *)
(*                (_Y_.)'  ._   )  .'_ ` ....-..-'         *)
(*               _..`--'_..-_/  /--'_.' ,'                 *)
(*              (il),-''  (li),'  ((!.-'                   *)
(*  ���������������������������������������������������ͻ  *)
(*  � TradeMark of ArashCorp Intl. & Real-Time systems. �  *)
(*  �     All rights reserved. Copyright(c) 1999.       �  *)
(*  ���������������������������������������������������ͼ  *)
(*                                                         *)
(*   Library Name   : FastGEO                              *)
(*   Production Date: 10/10/1999                           *)
(*   Library Version: 1.0.1                                *)
(*   Auther         : Arash Partow                         *)
(*   e-mail         : arashp@hotmail.com                   *)
(*                                                         *)
(* Capabilites:                                            *)
(*                                                         *)
(* 1. Fastest segment intersection detection routine,      *)
(* 2. Fastest point in triangle detection routine,         *)
(* 2. Fastest point in circle detection routine            *)
(* 4. Fastest Universal point in polygon detection routine *)
(* 5. Fastest Polygon segment intersection routine         *)
(* 6. Basic rectangular pair intersection routine          *)
(*                                                         *)
(* Areas which need further work:                          *)
(*                                                         *)
(* 1.) Need to add a test for point intersection or point  *)
(*     touching, for the PntInPolygon.                     *)
(*                                                         *)
(* 2.) Must add compatibility for non-ordered points of    *)
(*     a polygon.                                          *)
(*                                                         *)
(* 3.) Possibility of adding convex hull function, convex  *)
(*     polygons and simple types of concave polygons.      *)
(*                                                         *)
(* 4.) Add basic triangulation of convex polygons          *)
(*                                                         *)
(* 5.) Function to determine actual point of segments      *)
(*     intersection.                                       *)
(*                                                         *)
(* 6.) Collinear segment intersection test for non-complex *)
(*     segments                                            *)
(***********************************************************)

Unit FastGeo;

Interface


(* Basic geometric structure types *)

(* Vertex type *)
Type TPoint = Record x,y:LongInt; End;

(* Quadex type *)
Type TRect  = Array[1..2] Of TPoint;

(* Triangle type *)
Type TTriangle = Array[1..3] Of TPoint;

(* Segment type *)
Type TLine  = Array[1..2] Of TPoint;


Const
  (* Orientation constants *)
     RightHand        = -1;
     LeftHand         = +1;
     Collinear        = +0;
     Clockwise        = -1;
     CounterClockwise = +1;

(* Very fast and efficient 32-bit math code *)
Function LongAdd(A,B:LongInt):LongInt; Inline($66/$58/$66/$5B/$66/$03/$C3/$66/$0F/$A4/$C2/$10);
Function LongSub(A,B:LongInt):LongInt; Inline($66/$5B/$66/$58/$66/$2B/$C3/$66/$0F/$A4/$C2/$10);
Function LongMul(A,B:LongInt):LongInt; Inline($66/$58/$66/$5B/$66/$F7/$E3/$66/$0F/$A4/$C2/$10);
(* End of math functions *)

(* Core of all geometric operation *)
Function Orientation(x1,y1,x2,y2,x3,y3:LongInt):ShortInt;
Function SegmentIntersect(x1,y1,x2,y2,x3,y3,x4,y4:LongInt):Boolean;
Procedure GetBoundsRect(NumPoly:LongInt; Var Poly:Array Of TPoint; Var Rect:TRect);
(* End of core operations *)

(* Various geometric operations *)
Function RectangleIntersect(x1,y1,x2,y2,x3,y3,x4,y4:LongInt):Boolean;
Function Distance(x1,y1,x2,y2:LongInt):LongInt;
Function PntInTriangle(Px,Py,x1,y1,x2,y2,x3,y3:LongInt):Boolean;
Function PntInCircle(Px,Py,x,y,Radius:LongInt):Boolean;
Function PntInRectangle(Px,Py:LongInt; Rec:TRect):Boolean;
Function UNIPntInRectangle(Px,Py,x1,y1,x2,y2:LongInt):Boolean;
Function PntInQuadix(Px,Py,x1,y1,x2,y2,x3,y3,x4,y4:LongInt):Boolean;
(* End of geometric operations *)

(* Polygon based geometric functions *)
Function PolygonSegmentIntersect(x1,y1,x2,y2:LongInt; NumPoly:LongInt; Var Poly: Array Of TPoint):Boolean;
Function PolygonInPolygon(NumPoly1,NumPoly2:LongInt; Var Poly1,Poly2: Array Of TPoint):Boolean;
Function PolygonIntersect(NumPoly1,NumPoly2:LongInt; Var Poly1,Poly2: Array Of TPoint):Boolean;
Function PntInPolygon(Px,Py:LongInt; NumPoly:LongInt; Var Poly: Array Of TPoint):Boolean;
Function UNIPntInPolygon(Var Poly:Array of TPoint; Px,Py:LongInt; NumPoly:LongInt):Boolean;
(* End of polygon based algorithims *)

Implementation

(* Orientation procedures, is a basic vector product   *)
(* - value means orientation is to the right-hand side *)
(* + value means orientation is to the left-hand side  *)
(* 0 value menas both vectors are  identical           *)
Function Orientation(x1,y1,x2,y2,x3,y3:LongInt):ShortInt;
 Var Z: LongInt;
 Begin
  (* 32-Bit math code version of vector product    *)
  {Z:=LongAdd(LongAdd(LongMul(x1,LongSub(y2,y3)),LongMul(x2,LongSub(y3,y1))),LongMul(x3,LongSub(y1,y2)));}

  (* Nomral mathematical version of vector product *)
   Z:=(x1*(y2-y3))+(x2*(y3-y1))+(x3*(y1-y2));

  If Z < 0 Then Orientation:=-1; (* Orientaion is to the right-hand side *)
  If Z > 0 Then Orientation:=+1; (* Orientaion is to the left-hand side  *)
  If Z = 0 then Orientation:=+0; (* Orientaion is neutral                *)
 End;
(* End Orientation *)

(* Segment intersection function        *)
(* True if interesction occures         *)
(* False if intersection does not occur *)
Function SegmentIntersect(x1,y1,x2,y2,x3,y3,x4,y4:LongInt):Boolean;
 Begin
  SegmentIntersect:=False;
  If (Orientation(x1,y1,x2,y2,x3,y3) <> Orientation(x1,y1,x2,y2,x4,y4)) And
     (Orientation(x3,y3,x4,y4,x1,y1) <> Orientation(x3,y3,x4,y4,x2,y2)) Then SegmentIntersect:=True;
 End;
(* End SegmentIntersect *)

(* Point in triangle region function                *)
(* True if point is within triangle                 *)
(* False if point is not within triangle            *)
(* Explanation: This procedure traverses through    *)
(* all the points of the triangle, tesing for its   *)
(* orientaion towards the point in question, if     *)
(* all the orientations are the same then the point *)
(* is within the triangle otherwise point is not    *)
(* within the triangle.                             *)
Function PntInTriangle(Px,Py,x1,y1,x2,y2,x3,y3:LongInt):Boolean;
 Begin
  PntInTriangle:=False;
  If (Orientation(x1,y1,x2,y2,Px,Py) = Orientation(x2,y2,x3,y3,Px,Py)) And
     (Orientation(x2,y2,x3,y3,Px,Py) = Orientation(x3,y3,x1,y1,Px,Py)) Then PntInTriangle:=True;
 End;
(* End PntInTriangle *)

(* Distance between 2 points algorithim              *)
(* Note: The result of this function will alsways be *)
(* of integer type.                                  *)
Function Distance(x1,y1,x2,y2:LongInt):LongInt;
 Begin
  Distance:=Round(Sqrt(LongAdd(Sqr(x2-x1),Sqr(y2-y1))));
 End;
(* End Distance *)

(* Point within circular region function           *)
(* True if point is within circle                  *)
(* False if point is not within circle             *)
(* Explanation: The function finds the distance    *)
(* between the center of the circle and the point  *)
(* If the distance between the two is less than or *)
(* equal to the length of the radius the point is  *)
(* within the circle.                              *)
(* Further analysis of the function shows that if  *)
(* distance between the center and the point is    *)
(* equal to the radius then the point is on the    *)
(* circumfereance of the circle.                   *)
Function PntInCircle(Px,Py,x,y,Radius:LongInt):Boolean;
 Begin
  PntInCircle:=False;
  If Distance(Px,Py,x,y) <= Radius then PntInCircle:=True
 End;
(* End PntInCircle *)

(* Very basic point in rectangle region function   *)
(* True value means point is within rectangle      *)
(* False value means point is not within rectangle *)
Function UniPntInRectangle(Px,Py,x1,y1,x2,y2:LongInt):Boolean;
 Begin
  UniPntInRectangle:=False;
  (* Check for points existance with the x1,y1-x2,y2 values *)
  If (x1 <= Px) And (Px <= x2) And
     (y1 <= Py) And (Py <= y2) Then UniPntInRectangle:=True
 End;
(* End UniPntInRectangle *)

(* Point in rectangle procedure using the TRect type instead *)
(* of having individual coordinate assignments.              *)
Function PntInRectangle(Px,Py:LongInt; Rec:TRect):Boolean;
 Begin
  PntInRectangle:=False;
  If (Px >= Rec[1].x) and (Px <= Rec[2].x) And (Py >= Rec[1].y) and (Py <= Rec[2].y) Then PntInRectangle:=True
 End;
(* End PntInRectangle *)

(* Point in a quadix type of polygon function          *)
(* Where a quadix is a polygon with 4 vertices where   *)
(* all all sides may not be parallel to each other     *)
(* This function is faster than using the PntInPolygon *)
(* function.                                           *)
(* Note: The quadix must be a convex type of polygon   *)
Function PntInQuadix(Px,Py,x1,y1,x2,y2,x3,y3,x4,y4:LongInt):Boolean;
 Begin
  PntInQuadix:=False;
  If (Orientation(x1,y1,x2,y2,Px,Py) = Orientation(x2,y2,x3,y3,Px,Py)) And
     (Orientation(x2,y2,x3,y3,Px,Py) = Orientation(x3,y3,x4,y4,Px,Py)) And
     (Orientation(x3,y3,x4,y4,Px,Py) = Orientation(x4,y4,x1,y1,Px,Py)) Then PntInQuadix:=True;
 End;
(* End PntInQuadix *)

Function RectangleIntersect(x1,y1,x2,y2,x3,y3,x4,y4:LongInt):Boolean;
 Begin
  RectangleIntersect:=False;
  If ((x1 <= x4) And (x2 <= x3) And (y1 <= y4) And (y2 <= y3)) Then RectangleIntersect:=True;
 End;
(* End RectangleIntersect *)

Function PolygonSegmentIntersect(x1,y1,x2,y2:LongInt; NumPoly:LongInt; Var Poly: Array Of TPoint):Boolean;
 Var I:LongInt;
 Begin
  PolygonSegmentIntersect:=False;
  I:=NumPoly;
  If Not SegmentIntersect(Poly[I].x, Poly[I].y, Poly[1].x, Poly[1].y,x1,y1,x2,y2) then
  For I:= 1 to NumPoly-1 do If SegmentIntersect(Poly[I].x,Poly[I].y,Poly[I+1].x,Poly[I+1].y,x1,y1,x2,y2) Then
      Begin
       PolygonSegmentIntersect:=True;
       Break;
      End;
 End;
(* End PolygonSegmentIntersect *)

(* This function passes a true value of poly1 is within poly2 *)
(* or if poly2 is within poly1.                               *)
Function PolygonInPolygon(NumPoly1,NumPoly2:LongInt; Var Poly1,Poly2: Array Of TPoint):Boolean;
 Var I,Poly1Cnt,Poly2Cnt:LongInt;
 Begin
  Poly1Cnt:=0; Poly2Cnt:=0;
  PolygonInPolygon:=False;
  For I:= 1 to NumPoly1 do If PntInPolygon(Poly1[I].x,Poly1[I].y,NumPoly2,Poly2) Then Inc(Poly1Cnt);
  For I:= 1 to NumPoly2 do If PntInPolygon(Poly2[I].x,Poly2[I].y,NumPoly1,Poly1) Then Inc(Poly2Cnt);
  If (Poly1Cnt = NumPoly1) Or (Poly2Cnt = NumPoly2) Then PolygonInPolygon:=True;
 End;
(* End PolygonInPolygon *)


(* This procedures passes a true value if Poly1 intersects Poly2    *)
(* A few bugs are still being worked on. Should be finished by next *)
(* version.                                                         *)
Function PolygonIntersect(NumPoly1,NumPoly2:LongInt; Var Poly1,Poly2: Array Of TPoint):Boolean;
 Var I,J:LongInt;
 Begin
  PolygonIntersect:=False;
  For I:= 1 to NumPoly1-1 do
   For J:= 1 to NumPoly2-1 do
     If SegmentIntersect(Poly1[I].x,Poly1[I].y,Poly1[I+1].x,Poly1[I+1].y,
                         Poly2[J].x,Poly2[J].y,Poly2[J+1].x,Poly2[J+1].y) Then
        Begin
         PolygonIntersect:=True;
         Break;
        End;
 End;
(* End PolygonIntersect *)

(* This function is only to be used with convex polygons that have no     *)
(* indents, this function is suited very well for 3D-Engine type polygons *)
Function PntInPolygon(Px,Py:LongInt; NumPoly:LongInt; Var Poly: Array Of TPoint):Boolean;
 Var I:LongInt;
     Orin:ShortInt;
 Begin
  PntInPolygon:=True;
  Orin:= Orientation(Poly[NumPoly].x,Poly[NumPoly].y,Poly[1].x,Poly[1].y,Px,Py);
  For I:= 1 to NumPoly-1 do If Orin <> Orientation(Poly[I].x,Poly[I].y,Poly[I+1].x,Poly[I+1].y,Px,Py) Then
      Begin
       PntInPolygon:=False;
       Break;
      End;
 End;
(* End PntInPolygon *)

(* This function is the universal type of point in region test     *)
(* It is slower than the other type of point detection functions.  *)
(* It is advised that you triangulate your polygon and use the     *)
(* PntInTriangle function instead of this.                         *)
(* This function should only be used when trying to test for basic *)
(* types of point in polygon procedures.                           *)
Function UNIPntInPolygon(Var Poly:Array of TPoint; Px,Py:LongInt; NumPoly:LongInt):Boolean;
 Var FarRightPt : TPoint;
     BoundsRect : TRect;
     I, Intersects : LongInt;
 Begin
  UNIPntInPolygon:=False;
  GetBoundsRect(NumPoly,Poly,BoundsRect);
  If Not PntInRectangle(Px,Py,BoundsRect) Then Exit;
  FarRightPt.x := Px+1;
  FarRightPt.y := Py+BoundsRect[2].y+1;
  Intersects := 0;
  For I := 0 to NumPoly-2 do
  If SegmentIntersect(Poly[I].x, Poly[I].y,Poly[I+1].x,Poly[I+1].y,Px,Py,FarRightPt.x,FarRightPt.y) Then Inc(Intersects);
  I:=NumPoly-1;
  If SegmentIntersect(Poly[I].x, Poly[I].y, Poly[0].x, Poly[0].y,Px,Py,FarRightPt.x,FarRightPt.y) Then Inc(Intersects);
  If Intersects Mod 2 <> 0 Then UNIPntInPolygon:=True;
 End;
(* End UNIPntInPolygon *)

Function Maximum(A,B:LongInt):LongInt;
 Begin
  If A > B then Maximum := A
     Else Maximum := B;
 End;

Function Minimum(A,B:LongInt):LongInt;
 Begin
  If A <= B then Minimum := A
     Else Minimum := B;
 End;

Procedure GetBoundsRect(NumPoly:LongInt; Var Poly:Array Of TPoint; Var Rect:TRect);
 Var I,xMax,yMax,xMin,yMin:LongInt;
 Begin
  XMax := -MaxInt;
  YMax := -MaxInt;
  XMin := +MaxInt;
  YMin := +MaxInt;
  For I := 0 to NumPoly-1 do
     Begin
      xMax := Maximum(xMax,Poly[I].x);
      yMax := Maximum(yMax,Poly[I].y);
      xMin := Minimum(xMin,Poly[I].x);
      yMin := Minimum(yMin,Poly[I].y);
     End;
  Rect[1].x:=XMin; Rect[1].y:=YMin;
  Rect[2].x:=XMax; Rect[2].y:=YMax;
 End;


Begin
End.


Definitions:

 Convex polygon: Some of the functions above will only work porperly if
                 the polygon passed to the function is of a convex type
                 or in other words simple polygon.
                 To be a convex polygon the following must be seen:

                 1.) Any triangle is a convex polygon

                 2.) All the points must have a clockwise or
                     counter clockwise type of order

                 3.) A spline between any 3 consecutive points
                     must have its mean towards the outter or inner end of
                     the polygon, and not outside the range of the 2
                     furtherest points.

                 4.) None of the segments of the polygon may intersect
                     with each other

                 Functions that need to have a convex polygon passed will
                 be stated in the definition of the function.
                 All UNI type procedures or function do not need convex
                 polygons to be passed.


Note:
     The orienation method was taken from Carlos Moreno's 2 part
     article on 2D geometric operation in the 1998 November issue of
     C\C++ Users journal.
     Mr Moreno can be contacted: moreno@cyerglobe.net,
     or via the www at: http://www.mochima.com

     Other geometric implementations where made by Arash Partow. All
     Point in region operations have been tested for convex polygons
     But may work for some simple concave type polygons. Take note that
     the UniPntInPolygon will not work if any of the vertex segments of
     the polygon are intersecting, This will result in an indefinte answer.

     The math code was generate by Arash Partow and has been extensively
     test for all types of i086 type models for 32-bit numeric operation.



---------------------------------- DEMO --------------------------------------
                                 CUT HERE
Program FastGeo_Demo;

(************************************************************)
(*                                                          *)
(*       This code is governed by GPL                       *)
(*       You should NOT use this code                       *)
(*       if you are not familiar with                       *)
(*       the term "G.P.L"                                   *)
(*                                                          *)
(*              ("`-''-'").___..--''"`-._                   *)
(*                6_ 6  )   `-.  (    ).`-..,...'`)         *)
(*                (_Y_.)'  ._   )  .'_ ` ....-..-'          *)
(*               _..`--'_..-_/  /--'_.' ,'                  *)
(*              (il),-''  (li),'  ((!.-'                    *)
(*  ���������������������������������������������������ͻ   *)
(*  � TradeMark of ArashCorp Intl. & Real-Time systems. �   *)
(*  �     All rights reserved. Copyright(c) 1999.       �   *)
(*  ���������������������������������������������������ͼ   *)
(*  FastGeo v1.0 DEMO                                       *)
(*  Aim: To demonstrate some of the basic capabilites of    *)
(*       the FastGeo library.                               *)
(*       The demo constructs a near circular shaped polygon,*)
(*       then feeds a series of points on the screen, and   *)
(*       checks if they are with the polygon region.        *)
(************************************************************)

Uses Crt,Graph,FastGeo;

Const CirCoor = 6;

Var I,J:Integer;
    Circle : Array[0..CirCoor-1] Of TPoint;

Procedure RotateXY(Rangle:Integer; Ox,Oy:Integer; Var Nx,Ny:LongInt);
 Begin
  Nx:=Round((Ox*Cos(Rangle*Pi/180))-(Oy*Sin(Rangle*Pi/180)))+320;
  Ny:=Round((Ox*Sin(Rangle*Pi/180))+(Oy*Cos(Rangle*Pi/180)))+240;
 End;

Procedure MakeCircleCoors;
 Var Ang:Integer;
     I:Integer;
 Begin
  Ang:=0;
  For I:= 0 to CirCoor-1 do
    If Ang <= 360 then
    Begin
     RotateXY(180-Ang,0,100,Circle[I].x,Circle[I].y);
     Inc(Ang,360 Div CirCoor);
    End;
 End;

Procedure DrawCircle;
 Begin
  SetColor(10);
  For I:= 0 to CirCoor-2  do
    Line(Circle[I].x,Circle[I].y,Circle[I+1].x,Circle[I+1].y);
    Line(Circle[I+1].x,Circle[I+1].y,Circle[0].x,Circle[0].y);
 End;

Begin
 InitGraph(I,J,'D:\Bp\Bgi');
 MakeCircleCoors;
 DrawCircle;
 For I:= 220  to 420 do
 For J:= 140 to 340 do If UNIPntInPolygon(Circle,I,J,CirCoor) Then PutPixel(I,J,10);
 Readln;
End.

--------------7EC501A120D126B6D5D2597A--