forked from emorway-usgs/modflow6
-
Notifications
You must be signed in to change notification settings - Fork 0
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
feat(GeomUtil): add geometry module with point_in_polygon utility (MO…
- Loading branch information
Showing
9 changed files
with
198 additions
and
34 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,90 @@ | ||
| module TestGeomUtil | ||
| use KindModule, only: I4B, DP | ||
| use testdrive, only: error_type, unittest_type, new_unittest, check, test_failed | ||
| use GeomUtilModule, only: point_in_polygon | ||
| use ConstantsModule, only: LINELENGTH | ||
| implicit none | ||
| private | ||
| public :: collect_geomutil | ||
|
|
||
| contains | ||
|
|
||
| subroutine collect_geomutil(testsuite) | ||
| type(unittest_type), allocatable, intent(out) :: testsuite(:) | ||
| testsuite = [ & | ||
| new_unittest("point_in_polygon_rect", & | ||
| test_point_in_polygon_rect) & | ||
| ] | ||
| end subroutine collect_geomutil | ||
|
|
||
| subroutine test_point_in_polygon_rect(error) | ||
| type(error_type), allocatable, intent(out) :: error | ||
| real(DP), allocatable :: poly(:, :) | ||
|
|
||
| ! allocate and define polygon | ||
| allocate (poly(4, 2)) | ||
| ! vertices in clockwise order | ||
| poly(1, :) = (/0.0_DP, 0.0_DP/) | ||
| poly(2, :) = (/0.0_DP, 1.0_DP/) | ||
| poly(3, :) = (/1.0_DP, 1.0_DP/) | ||
| poly(4, :) = (/1.0_DP, 0.0_DP/) | ||
|
|
||
| ! points inside polygon | ||
| call check(error, point_in_polygon(0.99_DP, 0.01_DP, poly)) | ||
| call check(error, point_in_polygon(0.5_DP, 0.5_DP, poly)) | ||
| call check(error, point_in_polygon(0.0001_DP, 0.9999_DP, poly)) | ||
| if (allocated(error)) return | ||
|
|
||
| ! points outside polygon | ||
| call check(error, (.not. point_in_polygon(0.5_DP, 1.00001_DP, poly))) | ||
| call check(error, (.not. point_in_polygon(-0.5_DP, 34.0_DP, poly))) | ||
| if (allocated(error)) return | ||
|
|
||
| ! points on vertices | ||
| call check(error, point_in_polygon(0.0_DP, 0.0_DP, poly)) | ||
| call check(error, point_in_polygon(1.0_DP, 0.0_DP, poly)) | ||
| call check(error, point_in_polygon(0.0_DP, 1.0_DP, poly)) | ||
| call check(error, point_in_polygon(1.0_DP, 1.0_DP, poly)) | ||
| if (allocated(error)) return | ||
|
|
||
| ! points on faces | ||
| call check(error, point_in_polygon(0.0_DP, 0.5_DP, poly)) | ||
| call check(error, point_in_polygon(0.5_DP, 0.0_DP, poly)) | ||
| call check(error, point_in_polygon(1.0_DP, 0.5_DP, poly)) | ||
| call check(error, point_in_polygon(0.5_DP, 1.0_DP, poly)) | ||
| if (allocated(error)) return | ||
|
|
||
| ! vertices counter-clockwise | ||
| poly(1, :) = (/0.0_DP, 0.0_DP/) | ||
| poly(2, :) = (/1.0_DP, 0.0_DP/) | ||
| poly(3, :) = (/1.0_DP, 1.0_DP/) | ||
| poly(4, :) = (/0.0_DP, 1.0_DP/) | ||
|
|
||
| ! points inside polygon | ||
| call check(error, point_in_polygon(0.99_DP, 0.01_DP, poly)) | ||
| call check(error, point_in_polygon(0.5_DP, 0.5_DP, poly)) | ||
| call check(error, point_in_polygon(0.0001_DP, 0.9999_DP, poly)) | ||
| if (allocated(error)) return | ||
|
|
||
| ! points outside polygon | ||
| call check(error, (.not. point_in_polygon(0.5_DP, 1.00001_DP, poly))) | ||
| call check(error, (.not. point_in_polygon(-0.5_DP, 34.0_DP, poly))) | ||
| if (allocated(error)) return | ||
|
|
||
| ! points on vertices | ||
| call check(error, point_in_polygon(0.0_DP, 0.0_DP, poly)) | ||
| call check(error, point_in_polygon(1.0_DP, 0.0_DP, poly)) | ||
| call check(error, point_in_polygon(0.0_DP, 1.0_DP, poly)) | ||
| call check(error, point_in_polygon(1.0_DP, 1.0_DP, poly)) | ||
| if (allocated(error)) return | ||
|
|
||
| ! points on faces | ||
| call check(error, point_in_polygon(0.0_DP, 0.5_DP, poly)) | ||
| call check(error, point_in_polygon(0.5_DP, 0.0_DP, poly)) | ||
| call check(error, point_in_polygon(1.0_DP, 0.5_DP, poly)) | ||
| call check(error, point_in_polygon(0.5_DP, 1.0_DP, poly)) | ||
| if (allocated(error)) return | ||
|
|
||
| end subroutine test_point_in_polygon_rect | ||
|
|
||
| end module TestGeomUtil |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,65 @@ | ||
| module GeomUtilModule | ||
| use KindModule, only: I4B, DP | ||
| implicit none | ||
| private | ||
| public :: between, point_in_polygon | ||
| contains | ||
|
|
||
| !> @brief Check if a value is between two other values (inclusive). | ||
| logical function between(x, a, b) | ||
| real(DP), intent(in) :: x, a, b | ||
| between = (x >= a .and. x <= b .or. x <= a .and. x >= b) | ||
| end function between | ||
|
|
||
| !> @brief Check if a point is within a polygon. | ||
| !! Vertices and edge points are considered in. | ||
| !! Reference: https://stackoverflow.com/a/63436180/6514033 | ||
| logical function point_in_polygon(x, y, poly) | ||
| real(DP), intent(in) :: x | ||
| real(DP), intent(in) :: y | ||
| real(DP), allocatable, intent(in) :: poly(:, :) | ||
| integer(I4B) :: i, ii, num_verts | ||
| real(DP) :: xa, xb, ya, yb, c = 0.0_DP | ||
|
|
||
| point_in_polygon = .false. | ||
| num_verts = size(poly, 1) | ||
| xa = poly(1, 1) | ||
| ya = poly(1, 2) | ||
|
|
||
| do i = 0, num_verts + 1 | ||
| ii = mod(i, num_verts) + 1 | ||
| xb = poly(ii, 1) | ||
| yb = poly(ii, 2) | ||
|
|
||
| ! boundary cases | ||
| if ((x == xa .and. y == ya) .or. (x == xb .and. y == yb)) then | ||
| ! vertex point | ||
| point_in_polygon = .true. | ||
| else if (ya == yb .and. y == ya .and. between(x, xa, xb)) then | ||
| ! horizontal edge | ||
| point_in_polygon = .true. | ||
| ! if within vertical range, cast a ray | ||
| else if (between(y, ya, yb)) then | ||
| if (y == ya .and. yb >= ya .or. y == yb .and. ya >= yb) then | ||
| xa = xb | ||
| ya = yb | ||
| cycle | ||
| end if | ||
| ! cross product | ||
| c = (xa - x) * (yb - y) - (xb - x) * (ya - y) | ||
| ! boundary case | ||
| if (c == 0) then | ||
| point_in_polygon = .true. | ||
| ! standard intersection | ||
| else if ((ya < yb) .eqv. (c > 0)) then | ||
| point_in_polygon = .not. point_in_polygon | ||
| end if | ||
| end if | ||
|
|
||
| if (point_in_polygon) exit | ||
| xa = xb | ||
| ya = yb | ||
| end do | ||
| end function point_in_polygon | ||
|
|
||
| end module GeomUtilModule |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters