Find the intersection of two polylines

Problem/Question/Abstract:

How to find the intersection of two polylines

Answer:

Solve 1:

You have to intersect each polygon segment set which has a collision of their overlapping rectangles defined by the start and end point of each segment except neigboring segments. That means m-1 * n-1 segments are possible. To make a fast overlapping (collision) set I use a xy hash tree based on quadtree decomposition of the segments. Here is the code for the line intersection:

{XYIntersect Container Intersection for 2 dimensional segments
XII/2001 TriplexWare; Written by A.Weidauer

Abstract: Representation for 2-dimensional segment intersections
Author: Alexander Weidauer (alex.weidauer@huckfinn.de)
Created: December 2001
Lastmod: December 2001

The Unit delivers a 2-dimensional segment intersection for several objects
represented by basic data types for I/O}

unit UXYIntersect;

interface

uses
  UConst; {Basic datatype definitions. See further down the page.)

{The function checks a possible segmentation of two segments. The first is defined by
the coordinate set S1( P1(x1, y1): P2(x2, y2)) and the second is
defined S2( P3(x3, y3): P4(x4, y4)) where P1, P2, P3, P4 be the points.
OutX and OutY represent the intersection coordinates and are only valid if
the function turns back the value TRUE. If the segments are paralell the flag
is set to TRUE and if the the segmets are paralell and overlapping eachother
then OutX, OutY keeping the heavy point of the 4 coordinates sets.
In this case you have to check the intervall borders.
The solution of the intersection is NOT a point, it is a SEGMENT again.}

function Isec(x1, y1, x2, y2, x3, y3, x4, y4: TDouble; var OutX, OutY: TDouble;
  var ParallelFlag: TBoolean): TBoolean;

implementation

function Isec(x1, y1, x2, y2, x3, y3, x4, y4: TDouble; var OutX, OutY: TDouble;
  var ParallelFlag: TBoolean): TBoolean;
var
  delta, rmu, l1, l2, l3: TDouble;
begin
  OutY := 0;
  OutX := 0;
  ParallelFlag := False;
  x2 := x2 - x1;
  y2 := y2 - y1;
  x4 := x4 - x3;
  y4 := y4 - y3;
  delta := x2 * y4 - y2 * x4;
  {First case segments are paralell !}
  if abs(delta) < 1 E - 8 then
  begin
    ParallelFlag := False;
    x2 := x2 + x1;
    x4 := x4 + x3;
    y2 := y2 + y1;
    y4 := y4 + y3;
    l1 := sqrt((x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1));
    l2 := sqrt((x3 - x1) * (x3 - x1) + (y3 - y1) * (y3 - y1));
    l3 := sqrt((x3 - x2) * (x3 - x2) + (y3 - y2) * (y3 - y2));
    if (l1 = l2 - l3) or (l1 = l2 + l3) then
    begin
      Result := true;
      Parallelflag := true;
      OutX := (x1 + x2 + x3 + x4) / 4;
      OutY := (y1 + y2 + y3 + y4) / 4;
    end
    else
      Result := False;
    Exit;
  end;
  {End of parallel case}
  rmu := ((x3 - x1) * y4 - (y3 - y1) * x4) / delta;
  if (rmu > 1) or (rmu < 0) then
  begin
    isec := False;
    Exit;
  end;
  OutX := x1 + RMU * x2;
  Outy := y1 + RMU * y2;
  x2 := x2 + x1;
  x4 := x4 + x3;
  y2 := y2 + y1;
  y4 := y4 + y3;
  if x1 > x2 then
    SwapDouble(x1, x2);
  if y1 > y2 then
    SwapDouble(y1, y2);
  if x3 > x4 then
    SwapDouble(x3, x4);
  if y3 > y4 then
    SwapDouble(y3, y4);
  {Rangecheck of the solution}
  if (outx > x2) or (outx < x1) or (outx > x4) or (outx < x3) or (outy > y2)
    or (outy < y1) or (outy > y4) or (outy < y3) then
  begin
    Result := False;
    exit;
  end;
  Result := True;
end;

end.

{Basic Datatype and file extention definitions
XII/2001TriplaxWare; Written by A.Weidauer

Abstract: Basic Datatype and file extention definitions
Author: Alexander Weidauer (alex.weidauer@huckfinn.de)
Created: December 2001
Lastmod: December 2001

The Unit delivers basic data types for I/O and their file extentions.}

unit UConst;

interface

const
  {registry destination set}
  cStorageName = 'Software\tsWB';
  {Maximal read buffer size for blocked bufferd reading}
  cMaxReadBuffer = 32000;
  {Maximal write buffer size for blocked bufferd writing}
  cMaxWriteBuffer = 32000;

type
  {Type encapsulation Boolean}
  TBoolean = Boolean;
  {Type encapsulation String}
  TString = string;
  {Type encapsulation Byte 8 Bit}
  TByte = Byte;
  {Type encapsulation Word 16 Bit}
  TWord = Word;
  {Type encapsulation LongWord 32 Bit}
  TLongWord = LongWord;
  {Type encapsulation SmallInt signed 8 Bit}
  TInt08 = SmallInt;
  {Type encapsulation Integer signed 16 Bit}
  TInt16 = ShortInt;
  {Type encapsulation Integer signed 32 Bit}
  TInt32 = LongInt;
  {Type encapsulation Integer signed 32 Bit as common integer}
  TInteger = TInt32;
  {Type encapsulation Integer signed 64 Bit}
  TInt64 = Int64;
  {Type encapsulation Single 4 Byte}
  TSingle = Single;
  {Type encapsulation Real 6 Byte}
  TReal = Real48;
  {Type encapsulation Double 8 Byte}
  TDouble = Double;
  {Type encapsulation Double 10 Byte}
  TExtended = Extended;
  {Swap data if a > b for String}
procedure SwapString(var a, b: TString);
{Swap data if a > b for TByte}
procedure SwapByte(var a, b: TByte);
{Swap data if a > b for TWord}
procedure SwapWord(var a, b: TWord);
{Swap data if a > b for TLongWord}
procedure SwapLongWord(var a, b: TLongWord);
{Swap data if a > b for TInt08}
procedure SwapInt08(var a, b: TInt08);
{Swap data if a > b for TInt16}
procedure SwapInt16(var a, b: TInt16);
{Swap data if a > b for TInt32}
procedure SwapInt32(var a, b: TInt32);
{Swap data if a > b for TInt64}
procedure SwapInt64(var a, b: TInt64);
{Swap data if a > b for single}
procedure SwapSingle(var a, b: TSingle);
{Swap data if a > b for double}
procedure SwapDouble(var a, b: TDouble);
{Swap data if a > b for Extended}
procedure SwapExtended(var a, b: TExtended);

implementation

procedure SwapString(var a, b: TString);
var
  r: TString;
begin
  r := a;
  a := b;
  b := r;
end;

procedure SwapByte(var a, b: TByte);
var
  r: TByte;
begin
  r := a;
  a := b;
  b := r;
end;

procedure SwapWord(var a, b: TWord);
var
  r: TWord;
begin
  r := a;
  a := b;
  b := r;
end;

procedure SwapLongWord(var a, b: TLongWord);
var
  r: TLongWord;
begin
  r := a;
  a := b;
  b := r;
end;

procedure SwapInt08(var a, b: TInt08);
var
  r: TInt08;
begin
  r := a;
  a := b;
  b := r;
end;

procedure SwapInt16(var a, b: TInt16);
var
  r: TInt16;
begin
  r := a;
  a := b;
  b := r;
end;

procedure SwapInt32(var a, b: TInt32);
var
  r: TInt32;
begin
  r := a;
  a := b;
  b := r;
end;

procedure SwapInt64(var a, b: TInt64);
var
  r: TInt64;
begin
  r := a;
  a := b;
  b := r;
end;

procedure SwapSingle(var a, b: TSingle);
var
  r: TSingle;
begin
  r := a;
  a := b;
  b := r;
end;

procedure SwapDouble(var a, b: TDouble);
var
  r: TDouble;
begin
  r := a;
  a := b;
  b := r;
end;

procedure SwapExtended(var a, b: TExtended);
var
  r: TExtended;
begin
  r := a;
  a := b;
  b := r;
end;

end.


Solve 2:

This function will return the list of points found on a line from (x1,y1) to (x2,y2).
The procedure will calculate the points in the direction the line is drawn. For the line (x1,y1)------ --(x2,y2) or the line (x2,y2)-------(x1,y1) the first point in the list is always (x1,y1) and the last point in the list is always (x2, y2).
Points are calculated along the axis with the most change so that as many points as possible are created for the line.

// The point object
TPointFill = class
  X: Integer;
  Y: Integer;
end;

// ----------------------------------------------------------------------------
// GetLinePoints
// ----------------------------------------------------------------------------

function GetLinePoints(X1, Y1, X2, Y2: Integer): TList;
var
  ChangeInX, ChangeInY, i, MinX, MinY, MaxX, MaxY, LineLength: Integer;
  ChangingX: Boolean;
  Point: TPointFill;
  ReturnList, ReversedList: TList;
begin
  ReturnList := TList.Create;
  ReversedList := TList.Create;

  // Get the change in the X axis and the Max & Min X values
  if X1 > X2 then
  begin
    ChangeInX := X1 - X2;
    MaxX := X1;
    MinX := X2;
  end
  else
  begin
    ChangeInX := X2 - X1;
    MaxX := X2;
    MinX := X1;
  end;

  // Get the change in the Y axis and the Max & Min Y values
  if Y1 > Y2 then
  begin
    ChangeInY := Y1 - Y2;
    MaxY := Y1;
    MinY := Y2;
  end
  else
  begin
    ChangeInY := Y2 - Y1;
    MaxY := Y2;
    MinY := Y1;
  end;

  // Find out which axis has the greatest change
  if ChangeInX > ChangeInY then
  begin
    LineLength := ChangeInX;
    ChangingX := True;
  end
  else
  begin
    LineLength := ChangeInY;
    ChangingX := false;
  end;

  // If the x's match then the line changes only on the Y axis
  if X1 = X2 then
  begin
    // Loop thru the points on the list, lowest to highest.
    for i := MinY to MaxY do
    begin
      Point := TPointFill.Create;
      Point.X := X1;
      Point.Y := i;
      ReturnList.Add(Point);
    end;

    // If the point was started on the right and went to the left then
                reverse the list.
    if Y1 > Y2 then
    begin
      ReversedList := ReversePointOrder(ReturnList);
      ReturnList := ReversedList;
    end;
  end
    // If the x's match then the line changes only on the Y axis
  else if Y1 = Y2 then
  begin
    // Loop thru the points on the list, lowest to highest.
    for i := MinX to MaxX do
    begin
      Point := TPointFill.Create;
      Point.X := i;
      Point.Y := Y1;
      ReturnList.Add(Point);
    end;

    // If the point was started on the bottom and went to the top then reverse the list.
    if X1 > X2 then
    begin
      ReversedList := ReversePointOrder(ReturnList);
      ReturnList := ReversedList;
    end;
  end
    // The line is on an angle
  else
  begin
    // Add the first point to the list.
    Point := TPointFill.Create;
    Point.X := X1;
    Point.Y := Y1;
    ReturnList.Add(Point);

    // Loop thru the longest axis
    for i := 1 to (LineLength - 1) do
    begin
      Point := TPointFill.Create;
      // If we are moving on the x axis then get the related Y point.
      if ChangingX then
      begin
        Point.y := Round((ChangeInY * i) / ChangeInX);
        Point.x := i;
      end
        // otherwise we are moving on the y axis so get the related X point.
      else
      begin
        Point.y := i;
        Point.x := Round((ChangeInX * i) / ChangeInY);
      end;

      // if y1 is smaller than y2 then we are moving in a Top to Bottom direction.
      // we need to add y1 to get the next y value.
      if Y1 < Y2 then
        Point.y := Point.Y + Y1
          // otherwise we are moving in a Bottom to Top direction.
        // we need to subtract y1 to get the next y value.
      else
        Point.Y := Y1 - Point.Y;

      // if X1 is smaller than X2 then we are moving in a Left to Right direction
      // we need to add x1 to get the next x value
      if X1 < X2 then
        Point.X := Point.X + X1
          // otherwise we are moving in a Right to Left direction
        // we need to subtract x1 to get the next x value.
      else
        Point.X := X1 - Point.X;

      ReturnList.Add(Point);
    end;
    // Add the second point to the list.
    Point := TPointFill.Create;
    Point.X := X2;
    Point.Y := Y2;
    ReturnList.Add(Point);
  end;
  Result := ReturnList;
end;

// ----------------------------------------------------------------------------
// ReversePointOrder
// ----------------------------------------------------------------------------

function ReversePointOrder(LinePointList: TList): TList;
var
  i: integer;
  NewPointList: TList;
  CurrentPointFill: TPointFill;
begin
  NewPointList := TList.Create;
  i := LinePointList.Count - 1;

  while i > -1 do
  begin
    CurrentPointFill := TPointFill(LinePointList.Items[i]);
    NewPointList.Add(CurrentPointFill);
    dec(i);
  end;

  Result := NewPointList;
end;