Point Inside 3D Convex Polygon in Fortran
5.00/5 (2 votes)
An algorithm to determine if a point is inside a 3D convex polygon for a given polygon vertices in Fortran
Introduction
This is a Fortran version of the original C++ algorithm which can be found here.
Background
Please refer to the original C++ algorithm here.
Using the Code
The algorithm is wrapped into a Fortran DLL GeoProc. The main test program LASProc reads point cloud data from a LAS file, then writes all inside points into a new LAS file. A LAS file is able to be displayed with a freeware FugroViewer.
To consume the GeoProc module library, prepare a GeoPolygonProc object first like this:
polygon = GeoPolygon(verticesArray)
call procObj%InitGeoPolygonProc(polygon);
To check if a point is inside the polygon:
procObj%PointInside3DPolygon(x, y, z)
Here are all modules in DLL library GeoProc.dll:
Please refer to C++ Version for the module explanation and comments.
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 has no built-in support for List data structure to dynamically add new element to an unknown size array. Part of the reason is that Fortran is one of the HPC (High Performance Computing) languages, mainly used in math computing, i.e., solving differential equation with huge Array manipulation. The Array is supposed to be more memory efficient and access faster than List which is built on top of Array, the disadvantage of Array is its size or shape needs to be known before allocating memory to it, in some cases, a maximum Array size has to be calculated and utilized in Array declaration or allocation, if the problem requires a more flexible Array elements manipulation rather than Array element value computing, List data structure will be a better choice.
Here is the module Utility to implement a List data structure.
- Method
pushis to add an integer to 1d integer array. - Method
push2dis to add a fixed size 1d integer array to a 2d integer array. - Method
list2dContainsis to check if the 2d integer array contains the 1d integer array.
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
Here is the LAS process program to use the GeoProc DLL library:
When porting program from other language to Fortran, remember Fortran Start Index is 1, many of other languages use 0 as Start Index. This difference will require the changes for File IO, Array index, Loop counter, etc. Here is an example, the SEEK position is 96 in other language with 0 as Start Index, but Fortran uses 97 to start reading, for SEEK position record_loc, Fortran start reading from 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
Here is a module test program for all GeoProc DLL library modules.
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
Here is the output of the module test program:
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
Compile
Use GNU Fortran Compiler.
To compile DLL library 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
To compile test program:
gfortran -o test.exe test.f03 -I..\GeoProc\module -L..\GeoProc\lib -lGeoProc
To compile LAS process program:
gfortran -o LASProc.exe LASProc.f03 -I..\GeoProc\module -L..\GeoProc\lib -lGeoProc
History
- 10th February, 2016: Initial version date
References
- 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
History
- 10th February, 2016: Initial version