ugeometry.pas

Delphi的大数运算演示 pudn上大多是VC的 所以传个Delphi的

unit UGeometry;

interface

uses Windows;
type
  Tline=record
    p1,p2:TPoint; {starting and ending points of a line segment}
  end;
  PPResult=(PPoutside, PPInside, PPVertex, PPEdge, PPError);

  {integer distance bertween two points}
  function intdist(const p1,p2:TPoint):integer;

  {Make a line from two points}
  function Line(const p1,p2:TPoint):Tline;

  {Do two line segments intersect?}
  function  intersect(L1,L2:TLine; var pointonborder:boolean; var IP:TPoint):boolean;
  function Linesintersect(line1,line2:TLine):boolean;
  {Find the line from a given point which is perpendicular to a given line}
  {p1 of the returned line is the given point, p2 is the intersection point of the
   given line and the reutrned line}
  function PointPerpendicularLine(L:TLine; P:TPoint):TLine;

  {Return the perpendicular distance from a point to given line}
  function PerpDistance(L:TLine; P:TPoint):Integer;

  {Define the line through a given point intersecting a given line at a given angle}
  function AngledLineFromLine(L:TLine; P:TPoint; Dist:extended; alpha:extended):TLine;

  {Is a given point internal to a given polygon?}
  {Result is of type PPResult}
  {Possible result values are: PPoutside, PPInside, PPVertex, PPEdge, PPError}
  function PointInPoly(const p:TPoint; Points:array of TPoint):PPResult;

implementation

Uses math;

{********** Line ***********8}
function Line(const p1,p2:TPoint):Tline;
{Make a line from two points}
begin
  result.p1:=p1;
  result.p2:=p2;
end;

{************ IntDist *********}
function intdist(const p1,p2:TPoint):integer;
{Integer distance }
begin
  result:=round(sqrt((sqr(p1.x-p2.x)+sqr(p1.y-p2.y))));
end;


{******************* LinesIntersect ****************}
function Linesintersect(line1,line2:TLine):boolean;
   {local procedure: getequation}
   procedure getequation(line:TLine;var slope,intercept:extended);
   begin
      If line.p1.x<>line.p2.x
      then slope:=(line.p2.y-line.p1.y)/ (line.p2.x-line.p1.x)
      else slope:=1E100;
      intercept:=line.p1.y-slope*line.p1.x;
   end;

  function overlap(const x,y:extended; const line:TLine):boolean;
  {return true if passed x and y are within the range of the endpoints of the
   passed line}
  begin
    If (x>=min(line.p1.x,line.p1.x))
        and (x<=max(line.p1.x,line.p2.x))
        and (y>=min(line.p1.y,line.p2.y))
        and (y<=max(line.p1.y,line.p2.y))
    then result:=true
    else result:=false;
  end;


var
  m1,m2,b1,b2:extended;
  x,y:extended;

begin
  {Method -
     a. find the equations of the lines,
     b. find where they intersect
     c. if the point is between the line segment end points for both lines,
        then they do intersect, otherwise not.
  }
  result:=false;
  {general equation of line: y=mx+b}
  getequation(line1,m1,b1);
  getequation(line2,m2,b2);

  {intersection condition
     any point (x,y) on line1 satisfies y=m1*x+b1, the intersection
     point also satisfies the line2 equation y=m2*x+b2,
     so y=m1*x+b1=m2*x+b2

     solve for X -
       m1*x+b1=m2*x+b2
       x=(b2-b1)/(m1-m2)
  }
  if m1<>m2 then
  begin
    x:=round((b2-b1)/(m1-m2));
    if abs(m1) < abs(m2) then
    y:=round(m1*x+b1) {try to get y intercept from smallest slope}
    else y:=round(m2*x+b2); {but try the other equation}
    {for intersection,  x and y must lie between the endpoints of both lines}
    If      (x>=min(line1.p1.x,line1.p2.x))
        and (x<=max(line1.p1.x,line1.p2.x))
        and (x>=min(line2.p1.x,line2.p2.x))
        and (x<=max(line2.p1.x,line2.p2.x))
        and (y>=min(line1.p1.y,line1.p2.y))
        and (y<=max(line1.p1.y,line1.p2.y))
        and (y>=min(line2.p1.y,line2.p2.y))
        and (y<=max(line2.p1.y,line2.p2.y))
    then result:=true;
  end
  else if b1=b2 then  {slopes and intercepts are equal}
  begin  {lines are colinear }
    {if either end of line 1 is within the x and y range of line2, or
     either end of line2 is with the x,y range of line1, then
     then the lines overlap. For simplicity, we'll just call it an intersection}
    with line1.p1 do result:=overlap(x,y, line2);
    If not result then  with line1.p2 do result:=overlap(x,y,Line2);
    if not result then with line2.p1 do result:=overlap(x,y,line1);
    if not result then with line2.p2 do result:=overlap(x,y,line1);
  end;
  {otherwise, slopes are equal and intercepts unequal
     ==> parallel lines ==> no intersection}
end;


(*
{This version  of Intersect has a bug and sometimes reports incorrect results}
{Use above version until it is fixed}
*)
{*************** Intersect ***************}
function  intersect(L1,L2:TLine; var pointonborder:boolean; var IP:TPoint):boolean;
{Return true if line segments L1 and L2 intersect,
 also indicate if just touching and return the intersection point coordinates}

     {***************** SameSide ************}
     function sameside(L:TLine; p1,p2:TPoint;
            var pointonborder:boolean; var IP:TPoint):int64;
      {p1,p2 on same side of line  L =>result>0
       opposite sides => result <0
       a point on L => result=0 }
       {IP=point of intersection}
      var
        dx,dy,dx1,dy1,dx2,dy2:int64;
        t,b,r,s:double;
      begin
        dx:=L.p2.x-L.P1.x;
        dy:=L.p2.y-L.P1.y;
        dx1:=p1.x-L.p1.x;
        dy1:=p1.y-L.p1.y;
        dx2:=p2.x-L.p2.x;
        dy2:=p2.y-L.p2.y;
        result:=(dx*dy1-dy*dx1)*(dx*dy2-dy*dx2);
        IP.x:=0;
        IP.y:=0;

        t:=dx1*(p2.y-p1.y)-dy1*(p2.x-p1.x);

        b:=dx*(p2.y-p1.y)-dy*(p2.x-p1.x);

        if b<>0 then r:=t/b else r:=10;

        if (r>=0) and (r<=1) then

        begin

          t:=dx1*dy-dy1*dx;

          b:=dx*(p2.y-p1.y)-dy*(p2.x-p1.x);

          if b<>0 then s:=t/b else s:=10;

          if (s>=0) and (s<=1) then

          begin

            IP.x:=round(l.p1.x + r*dx);

            IP.y:=round(l.p1.y+r*dy);

          end;

        end;
        if ((dx<>0) or (dy<>0)) and (result=0) then pointonborder:=true
        else pointonborder:=false;
      end; {sameside}


  var
    a,b:int64;
    pb:boolean;
    IP1,IP2:TPoint;
  begin
    pointonborder:=false;
    a:=sameside(L1,L2.p1,L2.p2, pb, IP1);
    if pb then pointonborder:=true;
    b:=sameside(L2,L1.p1,L1.p2,pb, IP2);
    if pb then
    begin
      pointonborder:=true;
      result:=true;
    end
    else  result:=(a<0) and (b<0);
    IP:=Ip1;
  end;


{**************** PointPerpendicularLine **********}
function PointPerpendicularLine(L:TLine; P:TPoint):TLine;
{Define the line through point P and perpendicular to Line L}
var
  m1,m2,b1,b2:extended;
  rx:extended;
begin
  with L do   {get slope and intercept for Line L}
  begin
     if p1.x<>p2.x then {make sure slope is not infinite}
     m1:=(p2.y-p1.y)/(p2.x-p1.x)  else m1:=1e20;
     b1:=p1.y-m1*p1.x;
  end;
  with result do
  begin
    p1:=p;
    if m1<>0 then
    begin
      m2:=-1/m1; {slope of perpendicular line}
      b2:=p.y-m2*p.x; {intercept of perpendicular line}
      rx:=(b2-b1)/(m1-m2);
      p2.x:=round(rx); {intersection point of the two lines}
      p2.y:=round(m2*rx+b2);
    end
    else
    begin  {line 1 was horizontal so this one must be vertical}
      p2.x:=p1.x;
      p2.y:=L.p2.y;
    end;
  end;
end;

{************ PerpDistance ************}
function PerpDistance(L:TLine; P:TPoint):Integer;
{Define the line through point P and perpendicular to Line L}
var
  m1,m2,b1,b2:extended;
  rx:extended;
  line2:TLine;
begin
  with L do   {get slope and intercept for Line L}
  begin
     if p1.x<>p2.x then {make sure slope is not infinite}
     m1:=(p2.y-p1.y)/(p2.x-p1.x)  else m1:=1e20;
     b1:=p1.y-m1*p1.x;
  end;
  with line2 do
  begin
    p1:=p;
    if m1<>0 then
     begin
        m2:=-1/m1; {slope of perpendicular line}
        b2:=p.y-m2*p.x; {intercept of perpendicular line}
        rx:=(b2-b1)/(m1-m2);
        p2.x:=round(rx); {intersection point of the two lines}
        p2.y:=round(m2*rx+b2);
    end
    else
    begin  {line 1 was horizontal so this one must be vertical}
      p2.x:=p1.x;
      p2.y:=L.p2.y;
    end;
    if (p1.x=p2.x) and (p1.y=p2.y) then result:=0
    {make sure that p2 of line2 is actually on the line}
    else if not Linesintersect(L,Line2) then result:=1000000
    else result:=round(intdist(p1,p2));
  end;
end;

{**************** AngledLineFromLine **********}
 function AngledLineFromLine(L:TLine; P:TPoint; Dist:extended; alpha:extended):TLine;
{compute a line fron point, P, on line, L, for a specified distance, dist
 at angle, alpha.  }
var theta, newangle:extended;
begin
  with L do
  begin
     if p1.x<>p2.x then {make sure slope is not infinite}
     theta:=arctan2((p1.y-p2.y),(p2.x-p1.x))
     else {vertical line}
     if p2.y>p1.y then theta:=pi/2 else theta:=-pi/2;
  end;
  with result do
  begin
    p1:=p;
    newangle:=theta+alpha;
    p2.x:=p.x+round(dist*cos(newangle));
    p2.y:=p.y-round(dist*sin(newangle));
  end;
end;





{************* PointInPoly *************}
function PointInPoly(const p:TPoint; Points:array of Tpoint):PPResult;
{determine where point P lines in relation to the polygon defined by
 arrray "points". }

 var
  i,count:integer;
  lt,lp:TLine;
  IP:TPoint;
  ob:boolean;
  OK:Boolean;
  loopcount:integer;

      function between(p1,p2,p3:TPoint):boolean;
      {Is p1 within the rectangle formed by p2 and p3?}
      begin
        if     (p1.x<=max(p2.x,p3.x)) and (p1.x>=min(p2.x,p3.x))
           and (p1.y<=max(p2.y,p3.y)) and (p1.y>=min(p2.y,p3.y))
        then result:=true
        else result:=false;
      end;

  begin
    count:=0;

    result:=PPError;
    for i:=0 to high(Points){nbrpoints} do
    begin
      if (p.x=points[i].x) and (p.y=points[i].y) then
      begin
        result:=PPVertex;
        break;
      end;
    end;
    if result<> PPVertex then
    begin
      {Lt = extension of line from point to infinity" in the X direction}
      lt.p1:=p;
      lt.p2.y:=0;  lt.p2.x:=1000;
      loopcount:=0; {safety stop, to prevent infinite loops}
      repeat
        OK:=true;
        inc(loopcount);
        lt.p2.y:=lt.p2.y+10;  {we'll change the angle of the line from th point
                                 to infinity until, we get one that does not intersect
                                 any vertices of the polygon}
        //drawline(lt);
        for i:=0 to high(Points){nbrpoints} do
        begin
          {Set line Lp equal to each edge of the polygon}
          lp.p1:=points[i];
          If i<high(Points){nbrpoints} then lp.p2:=points[i+1] else lp.p2:=points[0];

          {See if the edge and the line intersect}
          if intersect(lt,lp,ob,IP) then
          begin
            if  not ob then inc(count)
            else
            begin  {could be confusion here - the point may really be on the border,
                or the polygon vertex could be on the extension of our line running
                east from the given point}
              {If our extension line goes througn a vertex, we'll change the
               direction until we find one that does not intersect a vertex}
              if (lp.p1.y=Ip.y) and (lp.p1.x=Ip.x) then
              begin
                ok:=false;
                break;
              end
              else if (Ip.x=p.x) and (ip.y=p.y) then
              begin  {point is neither inside nor outside - it's on an edge}
                result:=PPEdge;
                break;
              end;
            end;
          end;
        end;
      until OK or (loopcount>200);
          If loopcount>200 then result:=PPError
      else
      if result = PPError then
      begin
        if (count mod 2)=1 then result:=PPInside
        else result:=PPOutside;
      end;
    end;
  end;


end.