14,974,240 members
Articles / High Performance Computing
Article
Posted 10 Feb 2016

9.9K views

# Point Inside 3D Convex Polygon in Fortran

Rate me:
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:

FORTRAN
```polygon = GeoPolygon(verticesArray)
call procObj%InitGeoPolygonProc(polygon);```

To check if a point is inside the polygon:

FORTRAN
`procObj%PointInside3DPolygon(x, y, z)`

Here are all modules in DLL library GeoProc.dll:

GeoPoint.f03:

FORTRAN
```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

interface operator (+)
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

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 module GeoPoint_Module```

GeoVector.f03:

FORTRAN
```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

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:

FORTRAN
```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

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:

FORTRAN
```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:

FORTRAN
```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 build-in support for List data structure to dynamicly adding new element to an unknown size array. Part of reason is that Fortran is one of the HPC (High Performance Computing) languages, mainly is used in math computing, i.e. solving diffrencial equation with huge Array manipulation. The Array is supposed to be more memory efficient and access faster than List which is buit 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 push is to add an integer to 1d integer array.

Method push2d is to add a fixed size 1d integer array to a 2d integer array.

Method list2dContains is to check if the 2d integer array contains the 1d integer array.

Utility.f03:

FORTRAN
```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:

FORTRAN
```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,
! varable 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
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:

FORTRAN
```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 postion is 97

read(1, pos = 132) x_scale, y_scale, z_scale, x_offset, y_offset, z_offset

allocate(buf_record(record_len))

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(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 writting record count
write(2, pos = 108) insideCount

close(unit=1)
close(unit=2)
close(unit=3)

deallocate(buf_record)

end program LASProc```

Here is a module test program for all GeoProc DLL library modules.

test.f03:

FORTRAN
```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:

FORTRAN
```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

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 prgram:

`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

Feburary 10, 2016: Initial version date.

## Reference

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://www.cs.rpi.edu/~szymansk/OOF90/bugs.html

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