2020/08/17

判斷點是否在多邊形內的方法

給定一個由多個點的 list 產生的多邊形,判斷另一個點座標,是否有包含在該多邊形的圖形中。

方法是從給定點座標開始,往隨便一個方向射出一條射線(例如水平往右射線),看看穿過多少條邊。如果穿過偶數次,表示點在簡單多邊形外部;如果穿過奇數次,表示點在簡單多邊形內部。

不過,要另外處理,當射線穿過頂點、射線與邊重疊的狀況,也就是給定點座標,與某一條邊共線的狀況。

這兩個連結有對方法做更詳細的說明

Point in Polygon

How to check if a given point lies inside or outside a polygon?

另外有提供一個 Java 版的實作

// A Java program to check if a given point
// lies inside a given polygon
// Refer https://www.geeksforgeeks.org/check-if-two-given-line-segments-intersect/
// for explanation of functions onSegment(),
// orientation() and doIntersect()
class GFG
{

    // Define Infinite (Using INT_MAX
    // caused overflow problems)
    static int INF = 10000;

    static class Point
    {
        int x;
        int y;

        public Point(int x, int y)
        {
            this.x = x;
            this.y = y;
        }
    };

    // Given three colinear points p, q, r,
    // the function checks if point q lies
    // on line segment 'pr'
    static boolean onSegment(Point p, Point q, Point r)
    {
        if (q.x <= Math.max(p.x, r.x) &&
            q.x >= Math.min(p.x, r.x) &&
            q.y <= Math.max(p.y, r.y) &&
            q.y >= Math.min(p.y, r.y))
        {
            return true;
        }
        return false;
    }

    // To find orientation of ordered triplet (p, q, r).
    // The function returns following values
    // 0 --> p, q and r are colinear
    // 1 --> Clockwise
    // 2 --> Counterclockwise
    static int orientation(Point p, Point q, Point r)
    {
        int val = (q.y - p.y) * (r.x - q.x)
                - (q.x - p.x) * (r.y - q.y);

        if (val == 0)
        {
            return 0; // colinear
        }
        return (val > 0) ? 1 : 2; // clock or counterclock wise
    }

    // The function that returns true if
    // line segment 'p1q1' and 'p2q2' intersect.
    static boolean doIntersect(Point p1, Point q1,
                               Point p2, Point q2)
    {
        // Find the four orientations needed for
        // general and special cases
        int o1 = orientation(p1, q1, p2);
        int o2 = orientation(p1, q1, q2);
        int o3 = orientation(p2, q2, p1);
        int o4 = orientation(p2, q2, q1);

        // General case
        if (o1 != o2 && o3 != o4)
        {
            return true;
        }

        // Special Cases
        // p1, q1 and p2 are colinear and
        // p2 lies on segment p1q1
        if (o1 == 0 && onSegment(p1, p2, q1))
        {
            return true;
        }

        // p1, q1 and p2 are colinear and
        // q2 lies on segment p1q1
        if (o2 == 0 && onSegment(p1, q2, q1))
        {
            return true;
        }

        // p2, q2 and p1 are colinear and
        // p1 lies on segment p2q2
        if (o3 == 0 && onSegment(p2, p1, q2))
        {
            return true;
        }

        // p2, q2 and q1 are colinear and
        // q1 lies on segment p2q2
        if (o4 == 0 && onSegment(p2, q1, q2))
        {
            return true;
        }

        // Doesn't fall in any of the above cases
        return false;
    }

    // Returns true if the point p lies
    // inside the polygon[] with n vertices
    static boolean isInside(Point polygon[], int n, Point p)
    {
        // There must be at least 3 vertices in polygon[]
        if (n < 3)
        {
            return false;
        }

        // Create a point for line segment from p to infinite
        Point extreme = new Point(INF, p.y);

        // Count intersections of the above line
        // with sides of polygon
        int count = 0, i = 0;
        do
        {
            int next = (i + 1) % n;

            // Check if the line segment from 'p' to
            // 'extreme' intersects with the line
            // segment from 'polygon[i]' to 'polygon[next]'
            if (doIntersect(polygon[i], polygon[next], p, extreme))
            {
                // If the point 'p' is colinear with line
                // segment 'i-next', then check if it lies
                // on segment. If it lies, return true, otherwise false
                if (orientation(polygon[i], p, polygon[next]) == 0)
                {
                    return onSegment(polygon[i], p,
                                     polygon[next]);
                }

                count++;
            }
            i = next;
        } while (i != 0);

        // Return true if count is odd, false otherwise
        return (count % 2 == 1); // Same as (count%2 == 1)
    }

    // Driver Code
    public static void main(String[] args)
    {
        Point polygon1[] = {new Point(0, 0),
                            new Point(10, 0),
                            new Point(10, 10),
                            new Point(0, 10)};
        int n = polygon1.length;
        Point p = new Point(20, 20);
        if (isInside(polygon1, n, p))
        {
            System.out.println("Yes");
        }
        else
        {
            System.out.println("No");
        }
        p = new Point(5, 5);
        if (isInside(polygon1, n, p))
        {
            System.out.println("Yes");
        }
        else
        {
            System.out.println("No");
        }
        p = new Point(-1, 10);
        n = polygon1.length;
        if (isInside(polygon1, n, p))
        {
            System.out.println("Yes");
        }
        else
        {
            System.out.println("No");
        }


        Point polygon2[] = {new Point(0, 0),
            new Point(5, 5), new Point(5, 0)};
        p = new Point(3, 3);
        n = polygon2.length;
        if (isInside(polygon2, n, p))
        {
            System.out.println("Yes");
        }
        else
        {
            System.out.println("No");
        }
        p = new Point(5, 1);
        if (isInside(polygon2, n, p))
        {
            System.out.println("Yes");
        }
        else
        {
            System.out.println("No");
        }
        p = new Point(8, 1);
        if (isInside(polygon2, n, p))
        {
            System.out.println("Yes");
        }
        else
        {
            System.out.println("No");
        }
    }
}

以下的實作,是根據 Java 版本的內容,改成 erlang 版

-module(polygon).

%% API
-export([
  test/0,
  test2/0
]).

-record(point, {
  x :: float(),
  y :: float()
}).
-type(point() :: #point{}).

-define(INF, 10000.0).

%% 檢查 Q 是否在 線段 PR 上
%% check if point Q lise on line segment(PR)
-spec on_segment(P :: point(), Q :: point(), R :: point() ) -> boolean().
on_segment(P, Q, R) ->
  #point{x = Px, y = Py} = P,
  #point{x = Qx, y = Qy} = Q,
  #point{x = Rx, y = Ry} = R,

%%  io:format("P=~p, Q=~p, R=~p~n", [P, Q, R]),
%%
%%  io:format("Qx=~p, max(Px, Rx)=~p~n", [Qx, max(Px, Rx)]),
%%  io:format("Qx=~p, min(Px, Rx)=~p~n", [Qx, min(Px, Rx)]),
%%  io:format("Qy=~p, max(Py, Ry)=~p~n", [Qy, max(Py, Ry)]),
%%  io:format("Qy=~p, min(Py, Ry)=~p~n", [Qy, min(Py, Ry)]),
  case (Qx =< max(Px, Rx)) and (Qx >= min(Px, Rx)) and( Qy =< max(Py, Ry)) and (Qy >= min(Py, Ry)) of
    true ->
      true;
    _ ->
      false
  end.

%% 查詢 P, Q, R 的順序
%% 0: 三點共線, 1: clockwise 順時鐘, 2: counterclockwise 逆時鐘
-spec orientation(P :: point(), Q :: point(), R :: point() ) -> integer().
orientation(P, Q, R) ->
  #point{x = Px, y = Py} = P,
  #point{x = Qx, y = Qy} = Q,
  #point{x = Rx, y = Ry} = R,

  Val = (Qy - Py) * (Rx - Qx) - (Qx - Px) * (Ry - Qy),
  case Val == 0 of
    true ->
      0;
    false ->
      case Val > 0 of
        true ->
          1;
        _ ->
          2
      end
  end.

%% 確認 line(P1, Q1) 是否有跟 line(P2, Q2) 相交
-spec intersect(P1 :: point(), Q1 :: point(), P2 :: point(), Q2 :: point() ) -> boolean().
intersect(P1, Q1, P2, Q2) ->
  Orientation1 = orientation(P1, Q1, P2),
  Orientation2 = orientation(P1, Q1, Q2),
  Orientation3 = orientation(P2, Q2, P1),
  Orientation4 = orientation(P2, Q2, Q1),

  % general case
  case (Orientation1 /= Orientation2) and (Orientation3 /= Orientation4) of
    true ->
      true;
    _ ->
      %% Special Cases
      %% p1, q1 and p2 are colinear and p2 lies on segment p1q1
      case (Orientation1 == 0) and on_segment(P1, P2, Q1) of
        true ->
          true;
        _ ->
          %% p1, q1 and p2 are colinear and q2 lies on segment p1q1
          case (Orientation2 == 0) and on_segment(P1, Q2, Q1) of
            true ->
              true;
            _ ->
              %% p2, q2 and p1 are colinear and  p1 lies on segment p2q2
              case (Orientation3 == 0) and on_segment(P2, P1, Q2) of
                true ->
                  true;
                _ ->
                  %% p2, q2 and q1 are colinear and q1 lies on segment p2q2
                  case (Orientation4 == 0) and on_segment(P2, Q1, Q2) of
                    true ->
                      true;
                    _ ->
                      false
                  end
              end
          end
      end
  end.

-spec in_polygon(Polygon :: list(), P :: point() ) -> boolean().
in_polygon(Polygon, P) ->
  case length(Polygon) < 3 of
    true ->
      false;
    _ ->
      #point{x = _Px, y = Py} = P,
      %% 產生一個點,最後要跟 P 連成一條 Py 到 INF 的水平線段
      PInf = #point{x=?INF, y=Py},

      %% 計算 line(P, PInf) 跟所有多邊形的邊線的交點的數量
      %% Count intersections of the above line with sides of polygon

      % 分成第一個點, 跟其他點 兩個 list
      {PolygonHead, PolygonLast} = lists:split(1, Polygon),

%%      io:format("PolygonHead=~p, PolygonLast=~p~n", [PolygonHead, PolygonLast]),

      % 把 第一個點接到 PolygonLast 後面
      Polygon2 = lists:append(PolygonLast, PolygonHead),
      % 合併 Polygon, Polygon2 為新的 list, [{Polygon, Polygon2}]
      PolygonList = lists:zip(Polygon, Polygon2),

%%      io:format("PolygonList=~p~n", [PolygonList]),

      {CountRes, OnSegmentFlagRes, OnSegmentRes} = lists:foldl(fun({P1, P2}, {Count, OnSegmentFlag, OnSegment}) ->

%%        io:format("  lists:foldl P1=~p, P2=~p, intersect(P1, P2, P, PInf)=~p, orientation(P1, P, P2)=~p, on_segment(P1, P, P2)=~p~n", [P1, P2, intersect(P1, P2, P, PInf), orientation(P1, P, P2), on_segment(P1, P, P2)]),
        %% 判斷 (P1, P2), (P, PInf) 是否有交點
        case intersect(P1, P2, P, PInf) of
          true ->
            case orientation(P1, P, P2) == 0 of
              true ->
                %% 如果 P 跟 line(P1, P2) 共線,判斷是否 P 有在該線段上
                {Count, true, on_segment(P1, P, P2)};
              _ ->
                {Count + 1, OnSegmentFlag, OnSegment}
            end;
          _ ->
            {Count, OnSegmentFlag, OnSegment}
        end
                                                               end,
        {0, false, false}, PolygonList),

%%      io:format("CountRes=~p, OnSegmentFlagRes=~p, OnSegmentRes=~p, (CountRes rem 2)=~p~n", [CountRes, OnSegmentFlagRes, OnSegmentRes, CountRes rem 2]),
      %% 判斷交點數量是否為奇數
      %% Return true if count is odd, false otherwise
      case OnSegmentFlagRes of
        true ->
          OnSegmentRes;
        _ ->
          case (CountRes rem 2) == 1 of
            true ->
              true;
            _ ->
              false
          end
      end
  end.

test() ->
  P = #point{x = 1.0, y = 2.0},
  Q = #point{x = 2.0, y = 4.0},
  R = #point{x = 3.0, y = 6.0},
  S = #point{x = 4.0, y = 8.0},

  ResOnSegment = on_segment(P, Q, R),
  io:format("ResOnSegment=~p~n", [ResOnSegment]),

  ResOrientation = orientation(P, Q, R),
  io:format("ResOrientation=~p~n", [ResOrientation]),

  ResIntersect = intersect(P, Q, R, S),
  io:format("ResIntersect=~p~n", [ResIntersect]),
  ok.

test2() ->
  Polygon = [ #point{x = 0.0, y = 0.0}, #point{x = 10.0, y = 0.0}, #point{x = 10.0, y = 10.0}, #point{x = 0.0, y = 10.0} ],

  Res1 = in_polygon(Polygon, #point{x = 20.0, y = 20.0}),
  io:format("Res1=~p~n", [Res1]),

  Res2 = in_polygon(Polygon, #point{x = 5.0, y = 5.0}),
  io:format("Res2=~p~n", [Res2]),

  Res3 = in_polygon(Polygon, #point{x = -1.0, y = 10.0}),
  io:format("Res3=~p~n", [Res3]),

  %%%%%%%%%%%
  Polygon2 = [ #point{x = 0.0, y = 0.0}, #point{x = 5.0, y = 5.0}, #point{x = 5.0, y = 0.0} ],
  Res4 = in_polygon(Polygon2, #point{x = 3.0, y = 3.0}),
  io:format("Res4=~p~n", [Res4]),

  Res5 = in_polygon(Polygon2, #point{x = 5.0, y = 1.0}),
  io:format("Res5=~p~n", [Res5]),

  Res6 = in_polygon(Polygon2, #point{x = 8.0, y = 1.0}),
  io:format("Res6=~p~n", [Res6]),

  ok.

沒有留言:

張貼留言