Fortran 中点是否在 3D 凸多边形内
5.00/5 (2投票s)
一个算法，用于确定点是否在给定多边形顶点 Fortran 的 3D 凸多边形内。
引言
这是原始 C++ 算法的 Fortran 版本，可以在这里找到。
背景
请参考原始 C++ 算法 这里。
Using the Code
该算法被封装到一个 Fortran DLL GeoProc 中。主测试程序 LASProc 从 LAS 文件读取点云数据，然后将所有内部点写入一个新的 LAS 文件。LAS 文件可以使用免费软件 FugroViewer 显示。
要使用 GeoProc 模块库，首先准备一个 GeoPolygonProc 对象，像这样：
polygon = GeoPolygon(verticesArray)
call procObj%InitGeoPolygonProc(polygon);
要检查一个点是否在多边形内部：
procObj%PointInside3DPolygon(x, y, z)
以下是 DLL 库 GeoProc.dll 中的所有模块：
有关模块的解释和注释，请参考C++ 版本。
GeoPoint.f03:
module GeoPoint_Module
implicit none
! data member
type GeoPoint
real :: x
real :: y
real :: z
end type GeoPoint
! constructor
interface GeoPoint
module procedure new_GeoPoint
end interface
! operator overloading
interface operator (+)
procedure add
end interface operator (+)
contains
type(GeoPoint) function new_GeoPoint(x, y, z) result(pt)
implicit none
real, intent(in) :: x
real, intent(in) :: y
real, intent(in) :: z
pt%x = x
pt%y = y
pt%z = z
end function
type(GeoPoint) function add(p1, p2) result(pt)
implicit none
type(GeoPoint), intent(in) :: p1
type(GeoPoint), intent(in) :: p2
pt%x = p1%x + p2%x
pt%y = p1%y + p2%y
pt%z = p1%z + p2%z
end function add
end module GeoPoint_Module
GeoVector.f03:
module GeoVector_Module
use GeoPoint_Module
implicit none
! data member
type GeoVector
private
type(GeoPoint) :: p0
type(GeoPoint) :: p1
real :: x
real :: y
real :: z
end type GeoVector
! constructor
interface GeoVector
module procedure new_GeoVector
end interface
! operator overloading
interface operator (*)
procedure multiple
end interface operator (*)
contains
type(GeoVector) function new_GeoVector(p0, p1) result(vt)
implicit none
type(GeoPoint), intent(in) :: p0
type(GeoPoint), intent(in) :: p1
vt%p0 = p0
vt%p1 = p1
vt%x = p1%x - p0%x
vt%y = p1%y - p0%y
vt%z = p1%z - p0%z
end function
type(GeoVector) function multiple(v1, v2) result(vt)
implicit none
type(GeoVector), intent(in) :: v1
type(GeoVector), intent(in) :: v2
vt%x = v1%y * v2%z - v1%z * v2%y
vt%y = v1%z * v2%x - v1%x * v2%z
vt%z = v1%x * v2%y - v1%y * v2%x
vt%p0 = v1%p0
vt%p1 = vt%p0 + GeoPoint(vt%x, vt%y, vt%z);
end function multiple
real function get_x(vt) result(ret)
implicit none
type(GeoVector) :: vt
ret = vt%x
end function get_x
real function get_y(vt) result(ret)
implicit none
type(GeoVector) :: vt
ret = vt%y
end function get_y
real function get_z(vt) result(ret)
implicit none
type(GeoVector) :: vt
ret = vt%z
end function get_z
end module GeoVector_Module
GeoPlane.f03:
module GeoPlane_Module
use GeoPoint_Module
use GeoVector_Module
implicit none
! Plane Equation: a * x + b * y + c * z + d = 0
! data member
type GeoPlane
real :: a
real :: b
real :: c
real :: d
contains
procedure :: initGeoPlane
end type GeoPlane
! constructor
interface GeoPlane
module procedure new_GeoPlane
end interface
! operator overloading
interface operator (*)
procedure multiplePoint
end interface operator (*)
interface operator (-)
module procedure negative
end interface operator (-)
contains
subroutine initGeoPlane(this, p0, p1, p2)
implicit none
class(GeoPlane) :: this
type(GeoPoint) :: p0, p1, p2
type(GeoVector) :: u, v, n
v = GeoVector(p0, p1);
u = GeoVector(p0, p2);
n = u * v;
! normal vector
this%a = get_x(n);
this%b = get_y(n);
this%c = get_z(n);
this%d = -(this%a * p0%x + this%b * p0%y + this%c * p0%z);
end subroutine initGeoPlane
type(GeoPlane) function new_GeoPlane(a, b, c, d) result(pl)
implicit none
real, intent(in) :: a
real, intent(in) :: b
real, intent(in) :: c
real, intent(in) :: d
pl%a = a
pl%b = b
pl%c = c
pl%d = d
end function
real function multiplePoint(pl, pt) result(ret)
implicit none
type(GeoPlane), intent(in) :: pl
type(GeoPoint), intent(in) :: pt
ret = pt%x * pl%a + pt%y * pl%b + pt%z * pl%c + pl%d
end function multiplePoint
type(GeoPlane) function negative(this) result(pl)
implicit none
type(GeoPlane), intent(in) :: this
pl%a = -this%a
pl%b = -this%b
pl%c = -this%c
pl%d = -this%d
end function
end module GeoPlane_Module
GeoPolygon.f03:
module GeoPolygon_Module
use GeoPoint_Module
use Utility_Module
implicit none
! data member
type GeoPolygon
type(GeoPoint), dimension(:), allocatable :: pts
integer, dimension(:), allocatable :: idx
integer :: n
end type GeoPolygon
! constructor
interface GeoPolygon
module procedure new_GeoPolygon
end interface
contains
type(GeoPolygon) function new_GeoPolygon(ptsIn) result(this)
implicit none
type(GeoPoint), dimension(:), intent(in) :: ptsIn
integer :: i, isize
isize = size(ptsIn)
this%n = isize
allocate(this%idx(isize))
allocate(this%pts(isize))
do i = 1, isize
this%pts(i) = ptsIn(i)
this%idx(i) = i
end do
end function new_GeoPolygon
subroutine destructor(this)
implicit none
type(GeoPolygon) :: this
if (allocated(this % pts)) deallocate(this % pts)
if (allocated(this % idx)) deallocate(this % idx)
end subroutine
end module GeoPolygon_Module
GeoFace.f03:
module GeoFace_Module
use GeoPoint_Module
use Utility_Module
implicit none
! data member
type GeoFace
type(GeoPoint), dimension(:), allocatable :: pts
integer, dimension(:), allocatable :: idx
integer :: n
end type GeoFace
! constructor
interface GeoFace
module procedure new_GeoFace
end interface GeoFace
contains
type(GeoFace) function new_GeoFace(ptsIn, idxIn) result(this)
implicit none
type(GeoPoint), dimension(:), intent(in) :: ptsIn
integer, dimension(:), intent(in) :: idxIn
integer :: i, isize
isize = size(ptsIn)
this%n = isize
allocate(this%idx(isize))
allocate(this%pts(isize))
do i = 1, isize
this%pts(i) = ptsIn(i)
this%idx(i) = idxIn(i)
end do
end function new_GeoFace
subroutine destructor(this)
implicit none
type(GeoFace) :: this
if (allocated(this % pts)) deallocate(this % pts)
if (allocated(this % idx)) deallocate(this % idx)
end subroutine
end module GeoFace_Module
Fortran 没有内置对 List 数据结构的支持，无法动态地将新元素添加到未知大小的数组。部分原因是 Fortran 是 HPC（高性能计算）语言之一，主要用于数学计算，即使用大量的 Array 操作求解微分方程。Array 应该比基于 Array 构建的 List 更节省内存，访问速度更快，但 Array 的缺点是，在为其分配内存之前，需要知道它的大小或形状。在某些情况下，必须计算出最大的 Array 大小并在 Array 声明或分配中使用它。如果问题需要更灵活的 Array 元素操作，而不是 Array 元素值计算，那么 List 数据结构将是一个更好的选择。
这是实现 List 数据结构的 Utility 模块。
- 方法
push是将一个整数添加到一个一维整数数组。 - 方法
push2d是将一个固定大小的一维整数数组添加到一个二维整数数组。 - 方法
list2dContains是检查二维整数数组是否包含一维整数数组。
Utility.f03:
module Utility_Module
contains
! list: 1d array
! element: real number
subroutine push(list, element)
implicit none
integer :: i, isize
integer, intent(in) :: element
integer, dimension(:), allocatable, intent(inout) :: list
integer, dimension(:), allocatable :: clist
if(allocated(list)) then
isize = size(list)
allocate(clist(isize+1))
do i=1,isize
clist(i) = list(i)
end do
clist(isize+1) = element
deallocate(list)
call move_alloc(clist, list)
else
allocate(list(1))
list(1) = element
end if
end subroutine push
! list: a 2d array
! element: a 1d array
! all element must have same size
subroutine push2d(list, element)
implicit none
integer :: i, j, isize, esize
integer, dimension(:), intent(in) :: element
integer, dimension(:,:), allocatable, intent(inout) :: list
integer, dimension(:,:), allocatable :: clist
if(allocated(list)) then
esize = size(element)
isize = size(list)/esize;
allocate(clist(isize+1, esize))
do i=1,isize
do j=1, esize
clist(i,j) = list(i,j)
end do
end do
do i=1, esize
clist(isize+1, i) = element(i)
end do
deallocate(list)
call move_alloc(clist, list)
else
esize = size(element)
allocate(list(1, esize))
do i=1,esize
list(1, i) = element(i)
end do
end if
end subroutine push2d
subroutine sortArray(array)
implicit none
integer, dimension(:), intent(inout) :: array
integer :: i, j, isize
integer :: temp
isize = size(array)
if (isize .gt. 1) then
do i = 1, isize - 1
do j = i + 1, isize
if(array(i) > array(j)) then
temp = array(j)
array(j) = array(i)
array(i) = temp
end if
end do
end do
end if
end subroutine sortArray
! check if 2d array list contains first esize elements in 1d array element
! esize <= size(element)
logical function list2dContains(list, element, esize) result(isContains)
implicit none
integer :: i, j, isize
integer :: temp
integer, dimension(:), allocatable :: tempListPart, tempElement
integer :: esize
integer, dimension(:), intent(in) :: element
integer, dimension(:,:), allocatable, intent(in) :: list
isize = size(list)/esize
isContains = .false.
if ( (size(list) .ge. esize) .and. &
(esize .gt. 1) .and. &
(mod(size(list), esize) .eq. 0) ) then
allocate(tempListPart(esize))
allocate(tempElement(esize))
do i=1, esize
tempElement(i) = element(i)
end do
call sortArray(tempElement)
do i=1,isize
do j=1, esize
tempListPart(j) = list(i, j)
end do
call sortArray(tempListPart)
isContains = .true.
do j=1, esize
if (tempListPart(j) .ne. tempElement(j)) then
isContains = .false.
exit
end if
end do
if(isContains) exit
end do
deallocate(tempListPart)
deallocate(tempElement)
end if
end function list2dContains
end module Utility_Module
GeoPolygonProc.f03:
module GeoPolygonProc_Module
use GeoPoint_Module
use GeoVector_Module
use GeoPlane_Module
use GeoPolygon_Module
use GeoFace_Module
use Utility_Module
implicit none
real, parameter :: MaxUnitMeasureError = 0.001
! data member
type GeoPolygonProc
type(GeoPolygon) :: polygon
real :: x0
real :: x1
real :: y0
real :: y1
real :: z0
real :: z1
real :: maxError
type(GeoFace), dimension(:), allocatable :: Faces
type(GeoPlane), dimension(:), allocatable :: FacePlanes
integer :: NumberOfFaces
contains
procedure :: InitGeoPolygonProc
procedure :: SetBoundary
procedure :: SetMaxError
procedure :: SetFacePlanes
procedure :: PointInside3DPolygon
procedure :: UpdateMaxError
end type GeoPolygonProc
contains
subroutine InitGeoPolygonProc(this, polygon)
implicit none
class(GeoPolygonProc) :: this
type(GeoPolygon), intent(in) :: polygon
this%polygon = polygon
call SetBoundary(this)
call SetMaxError(this)
call SetFacePlanes(this)
end subroutine InitGeoPolygonProc
subroutine SetBoundary(this)
implicit none
class(GeoPolygonProc) :: this
integer :: i, n
n = this%polygon%n;
this%x0 = this%polygon%pts(1)%x
this%y0 = this%polygon%pts(1)%y
this%z0 = this%polygon%pts(1)%z
this%x1 = this%polygon%pts(1)%x
this%y1 = this%polygon%pts(1)%y
this%z1 = this%polygon%pts(1)%z
do i = 1, n
if (this%polygon%pts(i)%x < this%x0) then
this%x0 = this%polygon%pts(i)%x
end if
if (this%polygon%pts(i)%y < this%y0) then
this%y0 = this%polygon%pts(i)%y
end if
if (this%polygon%pts(i)%z < this%z0) then
this%z0 = this%polygon%pts(i)%z
end if
if (this%polygon%pts(i)%x > this%x1) then
this%x1 = this%polygon%pts(i)%x
end if
if (this%polygon%pts(i)%y > this%y1) then
this%y1 = this%polygon%pts(i)%y
end if
if (this%polygon%pts(i)%z > this%z1) then
this%z1 = this%polygon%pts(i)%z
end if
end do
end subroutine SetBoundary
subroutine SetMaxError(this)
implicit none
class(GeoPolygonProc) :: this
this%maxError = (abs(this%x0) + abs(this%x1) + &
abs(this%y0) + abs(this%y1) + &
abs(this%z0) + abs(this%z1)) / 6 * MaxUnitMeasureError;
end subroutine SetMaxError
subroutine SetFacePlanes(this)
implicit none
class(GeoPolygonProc) :: this
logical :: isNewFace
integer :: i, j, k, m, n, p, l, onLeftCount, onRightCount, &
numberOfFaces, maxFaceIndexCount
real :: dis
type(GeoPoint) :: p0, p1, p2, pt
type(GeoPlane) :: trianglePlane, facePlane
type(GeoFace) :: face
type(GeoPoint), dimension(:), allocatable :: pts
type(GeoPlane), dimension(:), allocatable :: fpOutward
integer, dimension(:), allocatable :: &
pointInSamePlaneIndex, verticeIndexInOneFace, faceVerticeIndexCount, idx
! vertices indexes 2d array for all faces,
! variable face index is major dimension, fixed total number of vertices as minor dimension
! vertices index is the original index value in the input polygon
integer, dimension(:,:), allocatable :: faceVerticeIndex
! face planes for all faces defined with outward normal vector
allocate(fpOutward(this%polygon%n))
! indexes of other points that are in same plane
! with the 3 face triangle plane point
allocate(pointInSamePlaneIndex(this%polygon%n - 3))
! vertice indexes in one face
maxFaceIndexCount = this%polygon%n -1
allocate(verticeIndexInOneFace(maxFaceIndexCount))
numberOfFaces = 0
do i = 1, this%polygon%n
! triangle point #1
p0 = this%polygon%pts(i)
do j = i + 1, this%polygon%n
! triangle point #2
p1 = this%polygon%pts(j)
do k = j + 1, this%polygon%n
! triangle point #3
p2 = this%polygon%pts(k)
call trianglePlane%initGeoPlane(p0, p1, p2)
onLeftCount = 0
onRightCount = 0
m = 0
do l = 1, this%polygon%n
! any point except the 3 triangle points
if (l .ne. i .and. l .ne. j .and. l .ne. k) then
pt = this%polygon%pts(l)
dis = trianglePlane * pt
! if point is in the triangle plane
! add it to pointInSamePlaneIndex
if (abs(dis) .lt. this%maxError ) then
m = m + 1
pointInSamePlaneIndex(m) = l
else
if (dis .lt. 0) then
onLeftCount = onLeftCount + 1
else
onRightCount = onRightCount + 1
end if
end if
end if
end do
n = 0
do p = 1, maxFaceIndexCount
verticeIndexInOneFace(p) = -1
end do
! This is a face for a CONVEX 3d polygon.
! For a CONCAVE 3d polygon, this maybe not a face.
if ((onLeftCount .eq. 0) .or. (onRightCount .eq. 0)) then
! add 3 triangle plane point index
n = n + 1
verticeIndexInOneFace(n) = i
n = n + 1
verticeIndexInOneFace(n) = j
n = n + 1
verticeIndexInOneFace(n) = k
! if there are other vertices in this triangle plane
! add them to the face plane
if (m .gt. 0) then
do p = 1, m
n = n + 1
verticeIndexInOneFace(n) = pointInSamePlaneIndex(p)
end do
end if
! if verticeIndexInOneFace is a new face,
! add it in the faceVerticeIndex list,
! add the trianglePlane in the face plane list fpOutward
!print *, n, size(verticeIndexInOneFace), size(faceVerticeIndex)
isNewFace = .not. list2dContains(faceVerticeIndex, &
verticeIndexInOneFace, maxFaceIndexCount)
if ( isNewFace ) then
numberOfFaces = numberOfFaces + 1
call push2d(faceVerticeIndex, verticeIndexInOneFace)
if (onRightCount .eq. 0) then
fpOutward(numberOfFaces) = trianglePlane
else if (onLeftCount .eq. 0) then
fpOutward(numberOfFaces) = -trianglePlane
end if
call push(faceVerticeIndexCount, n)
end if
else
! possible reasons:
! 1. the plane is not a face of a convex 3d polygon,
! it is a plane crossing the convex 3d polygon.
! 2. the plane is a face of a concave 3d polygon
end if
end do ! k loop
end do ! j loop
end do ! i loop
! Number of Faces
this%NumberOfFaces = numberOfFaces
allocate(this%FacePlanes(this%NumberOfFaces))
allocate(this%Faces(this%NumberOfFaces))
! loop faces
do i = 1, this%NumberOfFaces
! set FacePlanes
this%FacePlanes(i) = GeoPlane(fpOutward(i)%a, fpOutward(i)%b, &
fpOutward(i)%c, fpOutward(i)%d)
! actual vertices count in the face
n = faceVerticeIndexCount(i)
allocate(pts(n))
allocate(idx(n))
! loop face vertices
do j = 1, n
k = faceVerticeIndex(i, j)
pt = GeoPoint(this%polygon%pts(k)%x, this%polygon%pts(k)%y, this%polygon%pts(k)%z)
pts(j) = pt
idx(j) = k
end do
!set Faces
this%Faces(i) = GeoFace(pts, idx)
deallocate(pts)
deallocate(idx)
end do
deallocate(pointInSamePlaneIndex)
deallocate(verticeIndexInOneFace)
deallocate(faceVerticeIndex)
deallocate(faceVerticeIndexCount)
deallocate(fpOutward)
end subroutine SetFacePlanes
! main function to be called. check if a point is inside 3d polygon
logical function PointInside3DPolygon(this, x, y, z) result(ret)
implicit none
class(GeoPolygonProc) :: this
real, intent(in) :: x, y, z
integer i
real :: dis
! If the point is in the opposite half space with normal vector for all 6 faces,
! then it is inside of the 3D polygon
ret = .true.
do i = 1, this%NumberOfFaces
dis = this%FacePlanes(i) * GeoPoint(x, y, z)
! If the point is in the same half space with normal vector for any face of the 3D polygon,
! then it is outside of the 3D polygon
if (dis .gt. 0) then
ret = .false.
exit
end if
end do
end function PointInside3DPolygon
! update maxError attribute valu of GeoPolygonProc object
! maxError is used in SetFacePlanes to threshold a maximum distance to
! check the points nearby the triangle plane if being considered to be inside the triangle plane
! maxError default value is calculated from polygon boundary in SetMaxError
subroutine UpdateMaxError(this, maxError)
implicit none
class(GeoPolygonProc) :: this
real, intent(in) :: maxError
this%maxError = maxError
end subroutine UpdateMaxError
end module GeoPolygonProc_Module
这是使用 GeoProc DLL 库的 LAS 处理程序。
从其他语言移植程序到 Fortran 时，请记住 Fortran 的起始索引是 1，许多其他语言使用 0 作为起始索引。这种差异需要更改文件 IO、数组索引、循环计数器等。这里有一个例子，SEEK 位置在其他语言中是 96，使用 0 作为起始索引，但 Fortran 使用 97 开始读取，对于 SEEK 位置 record_loc，Fortran 从 record_loc + 1 开始读取。
read(1, pos = 97) offset_to_point_data_value
...
read(1, pos = record_loc + 1) xi, yi, zi
LASProc.f03:
program LASProc
use GeoPoint_Module
use GeoVector_Module
use GeoPlane_Module
use GeoPolygon_Module
use GeoFace_Module
use Utility_Module
use GeoPolygonProc_Module
implicit none
! variables of GeoPolygonProc
type(GeoPoint) :: p1, p2, p3, p4, p5, p6, p7, p8
type(GeoPoint), dimension(8) :: verticesArray
type(GeoPolygon) :: polygon
type(GeoPolygonProc) :: procObj
! variables of FILE IO
integer (kind=4) :: xi, yi, zi
integer (kind=4) :: offset_to_point_data_value, record_num
integer (kind=2) :: record_len
double precision :: x_scale, y_scale, z_scale, x_offset, y_offset, z_offset
character, dimension(:), allocatable :: buf_header, buf_record
! local variables
integer :: i, j, insideCount, record_loc
real :: x, y, z
! get GeoPolygonProc object
p1 = GeoPoint( - 27.28046, 37.11775, - 39.03485)
p2 = GeoPoint( - 44.40014, 38.50727, - 28.78860)
p3 = GeoPoint( - 49.63065, 20.24757, - 35.05160)
p4 = GeoPoint( - 32.51096, 18.85805, - 45.29785)
p5 = GeoPoint( - 23.59142, 10.81737, - 29.30445)
p6 = GeoPoint( - 18.36091, 29.07707, - 23.04144)
p7 = GeoPoint( - 35.48060, 30.46659, - 12.79519)
p8 = GeoPoint( - 40.71110, 12.20689, - 19.05819)
verticesArray(1) = p1
verticesArray(2) = p2
verticesArray(3) = p3
verticesArray(4) = p4
verticesArray(5) = p5
verticesArray(6) = p6
verticesArray(7) = p7
verticesArray(8) = p8
polygon = GeoPolygon(verticesArray)
call procObj%InitGeoPolygonProc(polygon);
! process LAS file
open(unit=1, file="..\..\LASInput\cube.las", status="OLD", access="STREAM")
open(unit=2, file="..\..\LASOutput\cube_inside.las", status="REPLACE", access="STREAM")
open(unit=3, file="..\..\LASOutput\cube_inside.txt", status="REPLACE", action = "write")
! Fortran Start Index is 1 while many of other languages Start Index is 0
! The code for Array, File IO, Loop, etc, need to change
! when porting code from other language to Fortran.
! Here is the example, in C/C++, the SEEK position is 96, but in Fortran,
! the read position is 97
read(1, pos = 97) offset_to_point_data_value
read(1, pos = 106) record_len
read(1, pos = 108) record_num
read(1, pos = 132) x_scale, y_scale, z_scale, x_offset, y_offset, z_offset
allocate(buf_header(offset_to_point_data_value))
allocate(buf_record(record_len))
! read LAS header
read(1, pos = 1) buf_header
! write LAS header
write(2) buf_header
insideCount = 0
do i = 0, record_num - 1
! seek position
record_loc = offset_to_point_data_value + record_len * i;
! Fortran Start Index is 1
! read position = seek position + 1
read(1, pos = record_loc + 1) xi, yi, zi
x = xi * x_scale + x_offset;
y = yi * y_scale + y_offset;
z = zi * z_scale + z_offset;
if (x > procObj%x0 .and. x < procObj%x1 .and. &
y > procObj%y0 .and. y < procObj%y1 .and. &
z > procObj%z0 .and. z < procObj%z1 ) then
if (procObj%PointInside3DPolygon(x, y, z)) then
! read LAS Point record
read(1, pos = record_loc + 1) buf_record
! write LAS Point record
write(2) buf_record
! write text LAS Point record
write(3, fmt='(F15.6, F15.6, F15.6)') x, y, z
insideCount = insideCount + 1
end if
end if
end do
! update record_len in LAS header with actual writing record count
write(2, pos = 108) insideCount
close(unit=1)
close(unit=2)
close(unit=3)
deallocate(buf_header)
deallocate(buf_record)
end program LASProc
这是一个所有 GeoProc DLL 库模块的模块测试程序。
test.f03:
program Test
use GeoPoint_Module
use GeoVector_Module
use GeoPlane_Module
use GeoPolygon_Module
use GeoFace_Module
use Utility_Module
use GeoPolygonProc_Module
implicit none
type(GeoPoint) :: p1, p2, p3, p4, p5, p6, p7, p8, pt, insidePoint, outsidePoint
type(GeoVector) :: vt, v1, v2
type(GeoPlane) :: pl
real :: dis
type(GeoPoint), dimension(8) :: verticesArray
type(GeoPolygon) :: polygon
integer :: i, j, n
type(GeoPoint), dimension(4) :: faceVerticesArray
integer, dimension(4) :: faceVerticeIndexArray
type(GeoFace) :: face
integer, dimension(:), allocatable :: arr1, arr2, arr3, arr4
integer, dimension(:,:), allocatable :: arr2d
logical :: b1, b2
type(GeoPolygonProc) :: procObj
p1 = GeoPoint( - 27.28046, 37.11775, - 39.03485)
p2 = GeoPoint( - 44.40014, 38.50727, - 28.78860)
p3 = GeoPoint( - 49.63065, 20.24757, - 35.05160)
p4 = GeoPoint( - 32.51096, 18.85805, - 45.29785)
p5 = GeoPoint( - 23.59142, 10.81737, - 29.30445)
p6 = GeoPoint( - 18.36091, 29.07707, - 23.04144)
p7 = GeoPoint( - 35.48060, 30.46659, - 12.79519)
p8 = GeoPoint( - 40.71110, 12.20689, - 19.05819)
print *, "GeoPoint Test:"
pt = p1 + p2
print *, pt%x, pt%y, pt%z
print *, "GeoVector Test:"
v1 = GeoVector(p1, p2)
v2 = GeoVector(p1, p3)
vt = v2 * v1
print *, get_x(vt), get_y(vt), get_z(vt)
print *, "GeoPlane Test:"
call pl%initGeoPlane(p1, p2, p3);
print *, pl%a, pl%b, pl%c, pl%d
pl = GeoPlane(1.0, 2.0, 3.0, 4.0);
print *, pl%a, pl%b, pl%c, pl%d
pl = -pl;
print *, pl%a, pl%b, pl%c, pl%d
dis = pl * GeoPoint(1.0, 2.0, 3.0);
print *, dis
print *, "GeoPolygon Test:"
verticesArray(1) = p1
verticesArray(2) = p2
verticesArray(3) = p3
verticesArray(4) = p4
verticesArray(5) = p5
verticesArray(6) = p6
verticesArray(7) = p7
verticesArray(8) = p8
polygon = GeoPolygon(verticesArray)
n = polygon%n
do i = 1, n
print *, polygon%idx(i), ": (", polygon%pts(i)%x, ", ", polygon%pts(i)%y, ", ", polygon%pts(i)%z, ")"
end do
print *, "GeoFace Test:"
faceVerticesArray(1) = p1
faceVerticesArray(2) = p2
faceVerticesArray(3) = p3
faceVerticesArray(4) = p4
faceVerticeIndexArray = (/ 1, 2, 3, 4 /);
face = GeoFace(faceVerticesArray, faceVerticeIndexArray);
n = face%n
do i = 1, n
print *, face%idx(i), ": (", face%pts(i)%x, ", ", face%pts(i)%y, ", ", face%pts(i)%z, ")"
end do
print *, "Utility Test:"
call push(arr1, 1)
call push(arr1, 2)
call push(arr1, 3)
call push(arr1, 4)
arr2 = (/ 4, 5, 6, 7 /)
arr3 = (/ 2, 3, 1, 4, 6 /)
arr4 = (/ 2, 3, 1, 5, 7 /)
call push2d(arr2d, arr1)
call push2d(arr2d, arr2)
b1 = list2dContains(arr2d, arr3, 4)
b2 = list2dContains(arr2d, arr4, 4)
print *, b1, b2
print *, "GeoPolygonProc Test:"
call procObj%InitGeoPolygonProc(polygon);
print *, procObj%x0, procObj%x1
print *, procObj%y0, procObj%y1
print *, procObj%z0, procObj%z1
print *, procObj%maxError
print *, procObj%NumberOfFaces
print *, "Face Planes:";
do i = 1, procObj%NumberOfFaces
print *, i, ": ", procObj%FacePlanes(i)%a, procObj%FacePlanes(i)%b, &
procObj%FacePlanes(i)%c, procObj%FacePlanes(i)%d
end do
print *, "Faces:"
do i = 1, procObj%NumberOfFaces
print *, "Face #", i, ":"
do j = 1, procObj%Faces(i)%n
print *, procObj%Faces(i)%idx(j), ": (", procObj%Faces(i)%pts(j)%x, &
procObj%Faces(i)%pts(j)%y, procObj%Faces(i)%pts(j)%z
end do
end do
insidePoint = GeoPoint(-28.411750, 25.794500, -37.969000)
outsidePoint = GeoPoint(-28.411750, 25.794500, -50.969000)
b1 = procObj%PointInside3DPolygon(insidePoint%x, insidePoint%y, insidePoint%z)
b2 = procObj%PointInside3DPolygon(outsidePoint%x, outsidePoint%y, outsidePoint%z)
print *, b1, ", ", b2
end program Test
这是模块测试程序的输出。
GeoPoint Test:
-71.6806030 75.6250153 -67.8234558
GeoVector Test:
-178.390884 160.813675 -319.868134
GeoPlane Test:
-178.390884 160.813675 -319.868134 -23321.6328
1.00000000 2.00000000 3.00000000 4.00000000
-1.00000000 -2.00000000 -3.00000000 -4.00000000
-18.0000000
GeoPolygon Test:
1 : ( -27.2804604 , 37.1177483 , -39.0348511 )
2 : ( -44.4001389 , 38.5072708 , -28.7886009 )
3 : ( -49.6306496 , 20.2475700 , -35.0516014 )
4 : ( -32.5109596 , 18.8580494 , -45.2978516 )
5 : ( -23.5914192 , 10.8173704 , -29.3044491 )
6 : ( -18.3609104 , 29.0770702 , -23.0414391 )
7 : ( -35.4805984 , 30.4665909 , -12.7951899 )
8 : ( -40.7111015 , 12.2068901 , -19.0581894 )
GeoFace Test:
1 : ( -27.2804604 , 37.1177483 , -39.0348511 )
2 : ( -44.4001389 , 38.5072708 , -28.7886009 )
3 : ( -49.6306496 , 20.2475700 , -35.0516014 )
4 : ( -32.5109596 , 18.8580494 , -45.2978516 )
Utility Test:
T F
GeoPolygonProc Test:
-49.6306496 -18.3609104
10.8173704 38.5072708
-45.2978516 -12.7951899
2.92348750E-02
6
Face Planes:
1 : -178.390884 160.813675 -319.868134 -23321.6328
2 : 104.610008 365.194000 125.259911 -5811.86768
3 : 342.393494 -27.7903938 -204.924881 2372.95679
4 : -342.393677 27.7906227 204.925003 -10372.9639
5 : -104.609970 -365.193939 -125.260048 -2188.13623
6 : 178.390854 -160.813873 319.868256 15321.6396
Faces:
Face # 1 :
1 : ( -27.2804604 37.1177483 -39.0348511
2 : ( -44.4001389 38.5072708 -28.7886009
3 : ( -49.6306496 20.2475700 -35.0516014
4 : ( -32.5109596 18.8580494 -45.2978516
Face # 2 :
1 : ( -27.2804604 37.1177483 -39.0348511
2 : ( -44.4001389 38.5072708 -28.7886009
6 : ( -18.3609104 29.0770702 -23.0414391
7 : ( -35.4805984 30.4665909 -12.7951899
Face # 3 :
1 : ( -27.2804604 37.1177483 -39.0348511
4 : ( -32.5109596 18.8580494 -45.2978516
5 : ( -23.5914192 10.8173704 -29.3044491
6 : ( -18.3609104 29.0770702 -23.0414391
Face # 4 :
2 : ( -44.4001389 38.5072708 -28.7886009
3 : ( -49.6306496 20.2475700 -35.0516014
7 : ( -35.4805984 30.4665909 -12.7951899
8 : ( -40.7111015 12.2068901 -19.0581894
Face # 5 :
3 : ( -49.6306496 20.2475700 -35.0516014
4 : ( -32.5109596 18.8580494 -45.2978516
5 : ( -23.5914192 10.8173704 -29.3044491
8 : ( -40.7111015 12.2068901 -19.0581894
Face # 6 :
5 : ( -23.5914192 10.8173704 -29.3044491
6 : ( -18.3609104 29.0770702 -23.0414391
7 : ( -35.4805984 30.4665909 -12.7951899
8 : ( -40.7111015 12.2068901 -19.0581894
T , F
编译
使用 GNU Fortran 编译器。
编译 DLL 库 GeoProc.dll
gfortran -c .\src\Utility.f03
gfortran -c .\src\GeoPoint.f03
gfortran -c .\src\GeoVector.f03 -I.
gfortran -c .\src\GeoPlane.f03 -I.
gfortran -c .\src\GeoPolygon.f03 -I.
gfortran -c .\src\GeoFace.f03 -I.
gfortran -c .\src\GeoPolygonProc.f03 -I.
gfortran -shared -mrtd -o GeoProc.dll *.o
编译测试程序
gfortran -o test.exe test.f03 -I..\GeoProc\module -L..\GeoProc\lib -lGeoProc
编译 LAS 处理程序
gfortran -o LASProc.exe LASProc.f03 -I..\GeoProc\module -L..\GeoProc\lib -lGeoProc
历史
- 2016年2月10日：初始版本日期
参考文献
- https://gcc.gnu.org/wiki/GFortran
- https://en.wikipedia.org/wiki/Fortran#Fortran_2003
- https://www.pgroup.com/lit/articles/insider/v3n1a3.htm
- http://fortranwiki.org/fortran/show/File+extensions
- http://stackoverflow.com/questions/28048508/how-to-add-new-element-to-dynamical-array-in-fortran90
- http://www.cs.rpi.edu/~szymansk/OOF90/bugs.html
- http://www.tek-tips.com/viewthread.cfm?qid=1572697
- http://stackoverflow.com/questions/10959386/fortran-array-of-variable-size-arrays
- https://gcc.gnu.org/onlinedocs/gfortran/SHAPE.html#SHAPE
- http://programmers.stackexchange.com/questions/221892/should-i-use-a-list-or-an-array
历史
- 2016年2月10日：初始版本