[空间几何] 如何实现判别一个点是否在四边形内

Program www_fcode_cn
  Implicit None
  Logical :: ins
  Call polygon_contains_point_2_2d(4, [0.,1.,1.,0.], 1.1, [1.,1.,0.,0.], 0.5, ins)
  Write (*, *) ins
End Program www_fcode_cn
Subroutine polygon_contains_point_2_2d(n, xn, xval, yn, yval, inside)
  ! *******************************************************************************
  ! ! POLYGON_CONTAINS_POINT_2_2D finds if a point is inside a convex polygon in 2D.
  ! Modified:
  ! 06 February 1999
  ! Author:
  ! John Burkardt
  ! Parameters:
  ! Input, integer N, the number of nodes or vertices in the polygon.
  ! N must be at least 3.
  ! Input, real XN(N), the X coordinates of the vertices.
  ! Input, real XVAL, the X coordinate of the point to be tested.
  ! Input, real YN(N), the Y coordinates of the vertices.
  ! Input, real YVAL, the Y coordinate of the point to be tested.
  ! Output, logical INSIDE, is .TRUE. if ( X,Y) is inside
  ! the polygon or on its boundary, and .FALSE. otherwise.
  Implicit None
  Integer n,i
  Logical inside
  Real x1,x2,x3
  Real xn(n),xval,y1,y2,y3,yn(n),yval
  inside = .False.
  ! A point is inside a convex polygon if and only if it is inside
  ! one of the triangles formed by X(1),Y(1) and any two consecutive
  ! points on the polygon's circumference.
  x1 = xn(1)
  y1 = yn(1)
  Do i = 2, n - 1
    x2 = xn(i)
    y2 = yn(i)
    x3 = xn(i+1)
    y3 = yn(i+1)
    Call triangle_contains_point_1_2d(x1, y1, x2, y2, x3, y3, xval, yval, inside)
    If (inside) Then
      Return
    End If
  End Do
  Return
End Subroutine polygon_contains_point_2_2d
Subroutine triangle_contains_point_1_2d(x1, y1, x2, y2, x3, y3, x, y, inside)
  ! *******************************************************************************
  ! ! TRIANGLE_CONTAINS_POINT_1_2D finds if a point is inside a triangle in 2D.
  ! Modified:
  ! 16 June 2001
  ! Author:
  ! John Burkardt
  ! Parameters:
  ! Input, real X1, Y1, X2, Y2, X3, Y3, the triangle vertices.
  ! The vertices should be given in counter clockwise order.
  ! Input, real X, Y, the point to be checked.
  ! Output, logical INSIDE, is .TRUE. if (X,Y) is inside
  ! the triangle or on its boundary, and .FALSE. otherwise.
  Implicit None
  Real c(3)
  Logical inside
  Real x,x1,x2,x3
  Real y,y1,y2,y3
  Call triangle_barycentric_2d(x1, y1, x2, y2, x3, y3, x, y, c)
  inside = .not.(any(c(1:3)<0.0E+00))
End Subroutine triangle_contains_point_1_2d
Subroutine triangle_barycentric_2d(x1, y1, x2, y2, x3, y3, x, y, c)
  ! *******************************************************************************
  ! ! TRIANGLE_BARYCENTRIC_2D finds the barycentric coordinates of a point in 2D.
  ! Discussion:
  ! The barycentric coordinate of point X related to vertex A can be
  ! interpreted as the ratio of the area of the triangle with
  ! vertex A replaced by vertex X to the area of the original
  ! triangle.
  ! Modified:
  ! 20 October 2001
  ! Author:
  ! John Burkardt
  ! Parameters:
  ! Input, real X1, Y1, X2, Y2, X3, Y3, the triangle vertices.
  ! The vertices should be given in counter clockwise order.
  ! Input, real X, Y, the point to be checked.
  ! Output, real C(3), the barycentric coordinates of (X,Y) with respect
  ! to the triangle.
  Implicit None
  Integer, Parameter :: n = 2
  Integer, Parameter :: nrhs = 1
  Real a(n, n+nrhs)
  Real c(3)
  Integer info
  Real x,x1,x2,x3
  Real y,y1,y2,y3
  ! Set up the linear system
  ! ( X2-X1  X3-X1 ) C1  = X-X1
  ! ( Y2-Y1  Y3-Y1 ) C2    Y-Y1
  ! which is satisfied by the barycentric coordinates of (X,Y).
  a(1, 1) = x2 - x1
  a(1, 2) = x3 - x1
  a(1, 3) = x - x1
  a(2, 1) = y2 - y1
  a(2, 2) = y3 - y1
  a(2, 3) = y - y1
  ! Solve the linear system.
  Call rmat_solve(a, n, nrhs, info)
  If (info/=0) Then
    Write (*, '(a)') ' '
    Write (*, '(a)') 'TRIANGLE_BARYCENTRIC_2D - Fatal error!'
    Write (*, '(a)') '  The linear system is singular.'
    Write (*, '(a)') '  The input data does not form a proper triangle.'
    Stop
  End If
  c(1) = a(1, 3)
  c(2) = a(2, 3)
  c(3) = 1.0E+00 - c(1) - c(2)
End Subroutine triangle_barycentric_2d
Subroutine rmat_solve(a, n, nrhs, info)
  ! *******************************************************************************
  ! ! RMAT_SOLVE uses Gauss-Jordan elimination to solve an N by N linear system.
  ! Modified:
  ! 08 November 2000
  ! Author:
  ! John Burkardt
  ! Parameters:
  ! Input/output, real A(N,N+NRHS), contains in rows and columns 1
  ! to N the coefficient matrix, and in columns N+1 through
  ! N+NRHS, the right hand sides.  On output, the coefficient matrix
  ! area has been destroyed, while the right hand sides have
  ! been overwritten with the corresponding solutions.
  ! Input, integer NRHS, the number of right hand sides.  NRHS
  ! must be at least 0.
  ! Output, integer INFO, singularity flag.
  ! 0, the matrix was not singular, the solutions were computed;
  ! J, factorization failed on step J, and the solutions could not
  ! be computed.
  Implicit None
  Integer n
  Integer nrhs
  Real a(n, n+nrhs)
  Real apivot,factor,temp
  Integer i,j,k
  Integer info,ipivot
  info = 0
  Do j = 1, n
    ipivot = j
    apivot = a(j, j)
    Do i = j + 1, n
      If (abs(a(i,j))>abs(apivot)) Then
        apivot = a(i, j)
        ipivot = i
      End If
    End Do
    If (apivot==0.0E+00) Then
      info = j
      Return
    End If
    Do i = 1, n + nrhs
      Call r_swap(a(ipivot,i), a(j,i))
    End Do
    a(j, j) = 1.0E+00
    a(j, j+1:n+nrhs) = a(j, j+1:n+nrhs)/apivot
    Do i = 1, n
      If (i/=j) Then
        factor = a(i, j)
        a(i, j) = 0.0E+00
        a(i, j+1:n+nrhs) = a(i, j+1:n+nrhs) - factor*a(j, j+1:n+nrhs)
      End If
    End Do
  End Do
End Subroutine rmat_solve
Subroutine r_swap(x, y)
  Implicit None
  Real x,y,z
  z = x
  x = y
  y = z
End Subroutine r_swap