Fortran Coder

查看: 10931|回复: 12
打印 上一主题 下一主题

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

[复制链接]

18

帖子

3

主题

0

精华

入门

F 币
85 元
贡献
48 点
跳转到指定楼层
楼主
发表于 2015-10-25 17:49:28 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式
如何实现使用Fortran判别一个点点是否在四边形内?能不能列举个简单的程序说明一下,比如叉乘判别法、面积判别法
分享到:  微信微信
收藏收藏 点赞点赞 点踩点踩

954

帖子

0

主题

0

精华

大师

F 币
184 元
贡献
75 点

规矩勋章元老勋章新人勋章水王勋章热心勋章

QQ
沙发
发表于 2015-10-25 19:58:21 | 只看该作者
本帖最后由 vvt 于 2015-10-25 20:00 编辑

John Burkardt 写过一个几何函数库,里面有很多类似的函数。2D的,3D的都有。
参考 http://fcode.cn/code_prof-11-1.html
我列几个函数给你,主要是调用 polygon_contains_point_2_2d ,其他函数是内部实用的。

这个函数可以判断一个二维的点,是否在一个(凸)多边形里。当多边形为四边形时,就是你想要的。

[Fortran] 纯文本查看 复制代码
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

18

帖子

3

主题

0

精华

入门

F 币
85 元
贡献
48 点
板凳
 楼主| 发表于 2015-10-26 18:52:53 | 只看该作者
本帖最后由 hipeilei 于 2015-10-31 15:28 编辑
vvt 发表于 2015-10-25 19:58
John Burkardt 写过一个几何函数库,里面有很多类似的函数。2D的,3D的都有。
参考 http://fcode.cn/code_p ...

非常感谢!后来我用面积法做出来了。

59

帖子

2

主题

0

精华

大师

F 币
810 元
贡献
476 点
地板
发表于 2015-11-9 23:43:07 | 只看该作者
四边形从对角线分成两个三角形,对于其中一个三角形,其内点与三个角点连线构成的三个小三角形面积之和等于这个三角的面积,如果其中一个小三角形面积等于零说明这个点在某条边上,如果三个小三角形的面积之和大于大三角形的面积,点在其外,已知三点坐标其面积有一个行列式公式计算其面积,编程应该容易

评分

参与人数 1贡献 +5 收起 理由
vvt + 5 很给力!

查看全部评分

136

帖子

3

主题

0

精华

版主

F 币
1964 元
贡献
1677 点

帅哥勋章管理勋章爱心勋章新人勋章热心勋章元老勋章

5#
发表于 2015-11-10 12:43:01 | 只看该作者
另外还有一种思路:以这个点链接四个端点,构成四个三角形,计算这四给三角形的面积和,与这个四边形面积做对比:若相等,则在其内或者边上,若大于四边形面积,则在四边形之外。

18

帖子

3

主题

0

精华

入门

F 币
85 元
贡献
48 点
6#
 楼主| 发表于 2015-11-12 10:11:57 | 只看该作者
kerb 发表于 2015-11-9 23:43
四边形从对角线分成两个三角形,对于其中一个三角形,其内点与三个角点连线构成的三个小三角形面积之和等于 ...

我就是这么做的,多谢

18

帖子

3

主题

0

精华

入门

F 币
85 元
贡献
48 点
7#
 楼主| 发表于 2015-11-12 10:12:33 | 只看该作者
aliouying 发表于 2015-11-10 12:43
另外还有一种思路:以这个点链接四个端点,构成四个三角形,计算这四给三角形的面积和,与这个四边形面积做 ...

我正是这么做的,我觉得这样的方法就已经满足我的要求了,多谢!

135

帖子

15

主题

0

精华

版主

F 币
1159 元
贡献
637 点

爱心勋章管理勋章

8#
发表于 2015-11-17 09:33:10 | 只看该作者
再取一个四边形内的点(如形心),现在有两个点,将这个两个点分别带入4条线的方程,如果两个点带入每条线的方程所得数的符号都一致,则点在多变形内。同样的方法可以判断,点在任意多变形内,或者多面体内。

136

帖子

3

主题

0

精华

版主

F 币
1964 元
贡献
1677 点

帅哥勋章管理勋章爱心勋章新人勋章热心勋章元老勋章

9#
发表于 2015-11-18 12:19:33 | 只看该作者
珊瑚虫 发表于 2015-11-17 09:33
再取一个四边形内的点(如形心),现在有两个点,将这个两个点分别带入4条线的方程,如果两个点带入每条线 ...

这个思想不错,不知道计算量怎么样,要不大家都写给小程序,一起测试测试?

135

帖子

15

主题

0

精华

版主

F 币
1159 元
贡献
637 点

爱心勋章管理勋章

10#
发表于 2015-11-18 15:20:51 | 只看该作者
aliouying 发表于 2015-11-18 12:19
这个思想不错,不知道计算量怎么样,要不大家都写给小程序,一起测试测试? ...

速度快,实现简单,精度高,判断点与多边形、多面体关系时的不二之选,值得拥有
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-4-25 08:16

Powered by Tencent X3.4

© 2013-2024 Tencent

快速回复 返回顶部 返回列表