Matlab中inpoly函数的fortran版本，支持向量（绝对原创，欢迎验证）

程序源码（发现问题及时联系yueliang2100@163.com谢谢）

module inpolygon
implicit none

contains
function inpoly(p,node,edge_,TOL_) result(in_on)

!  INPOLY: Point-in-polygon testing.
!
! Determine whether a series of points lie within the bounds of a polygon
! in the 2D plane. General non-convex, multiply-connected polygonal
! regions can be handled.
!
! SHORT SYNTAX:
!
!   in_on = inpoly(p,node);
!
!   p   : The points to be tested as an Nx2 array [x1 y1; x2 y2; etc].
!   node: The vertices of the polygon as an Mx2 array [X1 Y1; X2 Y2; etc].
!         The standard syntax assumes that the vertices are specified in
!         consecutive order.
!
!   in_on : An Nx1 logical array with IN_ON(:,1) = in, IN_ON(:,2) = on
!   in  : An Nx1 logical array with IN(i) = TRUE if P(i,:) lies within the
!         region.
!   on  : An Nx1 logical array with ON(i) = TRUE if P(i,:) lies on a
!        polygon edge. (A tolerance is used to deal with numerical
!        precision, so that points within a distance of
!        eps^0.8*norm(node(:),inf) from a polygon edge are considered "on"
!        the edge.
!
! LONG SYNTAX:
!
!  in_on = inpoly(p,node,edge);
!
!  edge: An Mx2 array of polygon edges, specified as connections between
!        the vertices in NODE: [n1 n2; n3 n4; etc]. The vertices in NODE
!        do not need to be specified in connsecutive order when using the
!        extended syntax.
!
! EXAMPLE:
!
!   PIPTest:       ! Will run a few examples
!
! The algorithm is based on the crossing number test, which counts the
! number of times a line that extends from each point past the right-most
! region of the polygon intersects with a polygon edge. Points with odd
! counts are inside. A simple implementation of this method requires each
! wall intersection be checked for each point, resulting in an O(N*M)
! operation count.
!
!   YUE HongLiang  : 2009-2010
!   Email          : yuehongliang0910@gmail.com
!   Last updated   : 11/02/2010 with CVF 6.6
!
! Problems or suggestions? Email me.
real,dimension(:,:),intent(in) :: p,node
integer,dimension(:,:),intent(in),optional :: edge_
real,intent(in),optional :: TOL_
logical,dimension(size(p,1),2) :: in_on
!--------------------------------------------------------------------------
logical,dimension(size(p,1)) :: in,on
integer::i,nnode,n,nc,n1,n2
real ::x1,x2,y1,y2,xmin,xmax,TOL,X(size(p,1)),Y(size(p,1)),delt(size(p,1))
integer,dimension(:,:),allocatable :: edge
!! ERROR CHECKING
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if ( present(TOL_) ) then
    TOL = TOL_
else
    TOL = 1.0e-12
end if
nnode = size(node,1)
if ( size(p,2)/=2 ) then
   write(*,'(A)') 'P must be an Nx2 array.'
   return
end if
if ( size(node,2)/=2 ) then
   write(*,'(A)') 'NODE must be an Mx2 array.'
   return
end if

if (present(edge_) ) then 
    allocate(edge(size(edge_,1),2))
    edge =edge_
    if ( size(edge,2)/=2 ) then
       write(*,'(A)') 'EDGE must be an Mx2 array.'
       return
    end if
    if ( (maxval(edge)>nnode).Or.(any(edge<1)) ) then
       write(*,'(A)') 'Invalid EDGE.' 
       return
    end if
else ! Build edge if not passed
    allocate(edge(nnode,2))
    edge = reshape( (/(i,i=1,nnode),(i,i=2,nnode),1/), (/nnode,2/) )
end if 

!! PRE-PROCESSING
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
in = .false.         ! Because we're dealing with mod(in,2) we don't have
                     ! to actually increment the crossing number, we can
                     ! just flip a logical at each intersection (faster!)
on = in
n  = size(p,1)
nc = size(edge,1)
! Choose the direction with the biggest range as the "y-coordinate" for the
! test. This should ensure that the sorting is done along the best
! direction for long and skinny problems wrt either the x or y axes.
X = p(:,1)
Y = p(:,2)
TOL = TOL*minval( maxval(p,dim=1) - minval(p,dim=1) )
!! MAIN LOOP
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Loop through edges
do i=1,nc 
    ! Nodes in current edge
    n1 = edge(i,1)
    n2 = edge(i,2)

    ! Endpoints - sorted so that [x1,y1] & [x2,y2] has y1<=y2
    !           - also get xmin = min(x1,x2), xmax = max(x1,x2)
    y1 = node(n1,2)
    y2 = node(n2,2)
    if (y1<y2) then
        x1 = node(n1,1)
        x2 = node(n2,1)
    else
        y1 = node(n2,2) ! y1 <=> y2
        y2 = node(n1,2)
        x1 = node(n2,1)
        x2 = node(n1,1);
    end if
    if (x1>x2) then
        xmin = x2
        xmax = x1
    else
        xmin = x1
        xmax = x2
    end if
    delt = abs( (Y-y1)*(X-x2)-(Y-y2)*(X-x1) )
    where( (Y>y1).and.(Y<=y2).and.((y2-y1)*(X-x1)<(Y-y1)*(x2-x1)) )             ! in the polygon 
!     where( (Y>=y1).and.(Y<y2).and.((y2-y1)*(X-x1)<(Y-y1)*(x2-x1)) ) 
        in = .not.in        
    elsewhere( (Y>=y1).and.(Y<=y2).and.(X>=xmin).and.(X<=xmax).and.(delt<TOL) ) ! on the edge
        on = on .or. .true. 
    end where
!     where( on.and.in ) in = .not.in
    in = (.not.(on.and.in).and.in)
end do
in_on(:,1) = in .or. on
in_on(:,2) = on
end function inpoly

end module inpolygon

测试代码

program PIPTest
use inpolygon
! real::node(3,2)=(/1.0,1.0,0.0,&
!                   0.0,1.0,1.0/),&
!       p(3,2)=(/1.0,0.6,2.0,&
!                0.5,0.6,1.0/)
! integer::edge(3,2)=(/1,2,3,&
!                      2,3,1/)
real::node(8,2)=(/1.0,0.0,-1.0, 0.0,2.0,0.0,-2.0, 0.0,&
                  0.0,1.0, 0.0,-1.0,0.0,2.0, 0.0,-2.0/),&
     p(16,2)=(/0.5,-0.5,-1.5,-0.5,1.5,-1.0,-3.0,-1.0,0.8,-0.8,-0.8, 0.8,0.4,-0.4,-0.4, 0.4,&
               0.0, 1.0, 0.0,-1.0,0.0, 2.0, 0.0,-2.0,0.8, 0.8,-0.8,-0.8,0.4, 0.4,-0.4,-0.4/)
integer::edge(8,2)=(/1,2,3,4,5,6,7,8,&
                     2,3,4,1,6,7,8,5/)
logical,dimension(16,2) :: inon
inon = inpoly(p,node,edge)
write(*,'(A5,16L)') 'in=',inon(:,1)
end program PIPTest