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
type GeoPoint
real :: x
real :: y
real :: z
end type GeoPoint
interface GeoPoint
module procedure new_GeoPoint
end interface
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
type GeoVector
private
type(GeoPoint) :: p0
type(GeoPoint) :: p1
real :: x
real :: y
real :: z
end type GeoVector
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:
module GeoPlane_Module
use GeoPoint_Module
use GeoVector_Module
implicit none
type GeoPlane
real :: a
real :: b
real :: c
real :: d
contains
procedure :: initGeoPlane
end type GeoPlane
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;
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
type GeoPolygon
type(GeoPoint), dimension(:), allocatable :: pts
integer, dimension(:), allocatable :: idx
integer :: n
end type GeoPolygon
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
type GeoFace
type(GeoPoint), dimension(:), allocatable :: pts
integer, dimension(:), allocatable :: idx
integer :: n
end type GeoFace
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
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:
module Utility_Module
contains
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
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
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
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
integer, dimension(:,:), allocatable :: faceVerticeIndex
allocate(fpOutward(this%polygon%n))
allocate(pointInSamePlaneIndex(this%polygon%n - 3))
maxFaceIndexCount = this%polygon%n -1
allocate(verticeIndexInOneFace(maxFaceIndexCount))
numberOfFaces = 0
do i = 1, this%polygon%n
p0 = this%polygon%pts(i)
do j = i + 1, this%polygon%n
p1 = this%polygon%pts(j)
do k = j + 1, this%polygon%n
p2 = this%polygon%pts(k)
call trianglePlane%initGeoPlane(p0, p1, p2)
onLeftCount = 0
onRightCount = 0
m = 0
do l = 1, this%polygon%n
if (l .ne. i .and. l .ne. j .and. l .ne. k) then
pt = this%polygon%pts(l)
dis = trianglePlane * pt
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
if ((onLeftCount .eq. 0) .or. (onRightCount .eq. 0)) then
n = n + 1
verticeIndexInOneFace(n) = i
n = n + 1
verticeIndexInOneFace(n) = j
n = n + 1
verticeIndexInOneFace(n) = k
if (m .gt. 0) then
do p = 1, m
n = n + 1
verticeIndexInOneFace(n) = pointInSamePlaneIndex(p)
end do
end if
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
end if
end do
end do
end do
this%NumberOfFaces = numberOfFaces
allocate(this%FacePlanes(this%NumberOfFaces))
allocate(this%Faces(this%NumberOfFaces))
do i = 1, this%NumberOfFaces
this%FacePlanes(i) = GeoPlane(fpOutward(i)%a, fpOutward(i)%b, &
fpOutward(i)%c, fpOutward(i)%d)
n = faceVerticeIndexCount(i)
allocate(pts(n))
allocate(idx(n))
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
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
logical function PointInside3DPolygon(this, x, y, z) result(ret)
implicit none
class(GeoPolygonProc) :: this
real, intent(in) :: x, y, z
integer i
real :: dis
ret = .true.
do i = 1, this%NumberOfFaces
dis = this%FacePlanes(i) * GeoPoint(x, y, z)
if (dis .gt. 0) then
ret = .false.
exit
end if
end do
end function PointInside3DPolygon
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
type(GeoPoint) :: p1, p2, p3, p4, p5, p6, p7, p8
type(GeoPoint), dimension(8) :: verticesArray
type(GeoPolygon) :: polygon
type(GeoPolygonProc) :: procObj
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
integer :: i, j, insideCount, record_loc
real :: x, y, z
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);
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")
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(1, pos = 1) buf_header
write(2) buf_header
insideCount = 0
do i = 0, record_num - 1
record_loc = offset_to_point_data_value + record_len * i;
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(2) buf_record
write(3, fmt='(F15.6, F15.6, F15.6)') x, y, z
insideCount = insideCount + 1
end if
end if
end do
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
History
- 10th February, 2016: Initial version
