Point Inside 3D Convex Polygon in Fortran

• Download source code - 5.3 MB

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 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
    
    ! 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

• 10^th 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

• 10^th February, 2016: Initial version