Skip to content

Commit e125e2f

Browse files
wpbonelliwpbonelli
authored
fix(GeomUtil): fix point_in_polygon, expand unit tests (MODFLOW-ORG#1423)
* loop over polygon edges was off by 2 * early termination logic was wrong
1 parent b305820 commit e125e2f

File tree

2 files changed

+215
-60
lines changed

2 files changed

+215
-60
lines changed

‎autotest/TestGeomUtil.f90

Lines changed: 203 additions & 48 deletions
Original file line numberDiff line numberDiff line change
@@ -1,90 +1,245 @@
11
module TestGeomUtil
22
use KindModule, only: I4B, DP
3-
use testdrive, only: error_type, unittest_type, new_unittest, check, test_failed
3+
use testdrive, only: check, error_type, new_unittest, test_failed, &
4+
to_string, unittest_type
45
use GeomUtilModule, only: point_in_polygon
56
use ConstantsModule, only: LINELENGTH
67
implicit none
78
private
89
public :: collect_geomutil
10+
private :: test_point_in_polygon
911

1012
contains
1113

1214
subroutine collect_geomutil(testsuite)
1315
type(unittest_type), allocatable, intent(out) :: testsuite(:)
1416
testsuite = [ &
15-
new_unittest("point_in_polygon_rect", &
16-
test_point_in_polygon_rect) &
17+
new_unittest("point_in_polygon_sq", &
18+
test_point_in_polygon_sq), &
19+
new_unittest("point_in_polygon_tri", &
20+
test_point_in_polygon_tri), &
21+
new_unittest("point_in_polygon_irr", &
22+
test_point_in_polygon_irr) &
1723
]
1824
end subroutine collect_geomutil
1925

20-
subroutine test_point_in_polygon_rect(error)
26+
subroutine test_point_in_polygon(error, shape, &
27+
poly, in_pts, out_pts, vert_pts, face_pts)
28+
type(error_type), allocatable, intent(inout) :: error
29+
character(len=*), intent(in) :: shape
30+
real(DP), allocatable, intent(in) :: poly(:, :)
31+
real(DP), allocatable, intent(in) :: in_pts(:, :)
32+
real(DP), allocatable, intent(in) :: out_pts(:, :)
33+
real(DP), allocatable, intent(in) :: vert_pts(:, :)
34+
real(DP), allocatable, intent(in) :: face_pts(:, :)
35+
integer(I4B) :: i
36+
real(DP) :: x, y
37+
38+
! test inside points
39+
do i = 1, size(in_pts, 1)
40+
x = in_pts(i, 1)
41+
y = in_pts(i, 2)
42+
call check(error, point_in_polygon(x, y, poly), &
43+
"point inside "//shape//" failed: " &
44+
//to_string(x)//", "//to_string(y))
45+
if (allocated(error)) return
46+
end do
47+
48+
! test outside points
49+
do i = 1, size(out_pts, 1)
50+
x = out_pts(i, 1)
51+
y = out_pts(i, 2)
52+
call check(error, (.not. point_in_polygon(x, y, poly)), &
53+
"point outside "//shape//" failed: " &
54+
//to_string(x)//", "//to_string(y))
55+
if (allocated(error)) return
56+
end do
57+
58+
! test vertex points
59+
do i = 1, size(vert_pts, 1)
60+
x = vert_pts(i, 1)
61+
y = vert_pts(i, 2)
62+
call check(error, point_in_polygon(x, y, poly), &
63+
"point on "//shape//" vertex failed: " &
64+
//to_string(x)//", "//to_string(y))
65+
if (allocated(error)) return
66+
end do
67+
68+
! test face points
69+
do i = 1, size(face_pts, 1)
70+
x = face_pts(i, 1)
71+
y = face_pts(i, 2)
72+
call check(error, point_in_polygon(x, y, poly), &
73+
"point on "//shape//" face failed: " &
74+
//to_string(x)//", "//to_string(y))
75+
if (allocated(error)) return
76+
end do
77+
end subroutine test_point_in_polygon
78+
79+
!> @brief Test a unit square
80+
subroutine test_point_in_polygon_sq(error)
2181
type(error_type), allocatable, intent(out) :: error
2282
real(DP), allocatable :: poly(:, :)
83+
real(DP), allocatable :: in_pts(:, :)
84+
real(DP), allocatable :: out_pts(:, :)
85+
real(DP), allocatable :: vert_pts(:, :)
86+
real(DP), allocatable :: face_pts(:, :)
2387

24-
! allocate and define polygon
2588
allocate (poly(4, 2))
26-
! vertices in clockwise order
89+
90+
allocate (in_pts(3, 2))
91+
in_pts(1, :) = (/0.99_DP, 0.01_DP/)
92+
in_pts(2, :) = (/0.5_DP, 0.5_DP/)
93+
in_pts(3, :) = (/0.0001_DP, 0.9999_DP/)
94+
95+
allocate (out_pts(2, 2))
96+
out_pts(1, :) = (/0.5_DP, 1.00001_DP/)
97+
out_pts(2, :) = (/-0.5_DP, 34.0_DP/)
98+
99+
allocate (vert_pts(4, 2))
100+
vert_pts(1, :) = (/0.0_DP, 0.0_DP/)
101+
vert_pts(2, :) = (/1.0_DP, 0.0_DP/)
102+
vert_pts(3, :) = (/0.0_DP, 1.0_DP/)
103+
vert_pts(4, :) = (/1.0_DP, 1.0_DP/)
104+
105+
allocate (face_pts(4, 2))
106+
face_pts(1, :) = (/0.0_DP, 0.5_DP/)
107+
face_pts(2, :) = (/0.5_DP, 0.0_DP/)
108+
face_pts(3, :) = (/1.0_DP, 0.5_DP/)
109+
face_pts(4, :) = (/0.5_DP, 1.0_DP/)
110+
27111
poly(1, :) = (/0.0_DP, 0.0_DP/)
28112
poly(2, :) = (/0.0_DP, 1.0_DP/)
29113
poly(3, :) = (/1.0_DP, 1.0_DP/)
30114
poly(4, :) = (/1.0_DP, 0.0_DP/)
31-
32-
! points inside polygon
33-
call check(error, point_in_polygon(0.99_DP, 0.01_DP, poly))
34-
call check(error, point_in_polygon(0.5_DP, 0.5_DP, poly))
35-
call check(error, point_in_polygon(0.0001_DP, 0.9999_DP, poly))
115+
call test_point_in_polygon(error, "clockwise square", &
116+
poly, in_pts, out_pts, vert_pts, face_pts)
36117
if (allocated(error)) return
37118

38-
! points outside polygon
39-
call check(error, (.not. point_in_polygon(0.5_DP, 1.00001_DP, poly)))
40-
call check(error, (.not. point_in_polygon(-0.5_DP, 34.0_DP, poly)))
119+
poly(1, :) = (/0.0_DP, 0.0_DP/)
120+
poly(2, :) = (/1.0_DP, 0.0_DP/)
121+
poly(3, :) = (/1.0_DP, 1.0_DP/)
122+
poly(4, :) = (/0.0_DP, 1.0_DP/)
123+
call test_point_in_polygon(error, "counter-clockwise square", &
124+
poly, in_pts, out_pts, vert_pts, face_pts)
41125
if (allocated(error)) return
42126

43-
! points on vertices
44-
call check(error, point_in_polygon(0.0_DP, 0.0_DP, poly))
45-
call check(error, point_in_polygon(1.0_DP, 0.0_DP, poly))
46-
call check(error, point_in_polygon(0.0_DP, 1.0_DP, poly))
47-
call check(error, point_in_polygon(1.0_DP, 1.0_DP, poly))
48-
if (allocated(error)) return
127+
deallocate (poly)
128+
deallocate (in_pts)
129+
deallocate (out_pts)
130+
deallocate (vert_pts)
131+
deallocate (face_pts)
132+
end subroutine test_point_in_polygon_sq
133+
134+
!> @brief Test a right triangle
135+
subroutine test_point_in_polygon_tri(error)
136+
type(error_type), allocatable, intent(out) :: error
137+
real(DP), allocatable :: poly(:, :)
138+
real(DP), allocatable :: in_pts(:, :)
139+
real(DP), allocatable :: out_pts(:, :)
140+
real(DP), allocatable :: vert_pts(:, :)
141+
real(DP), allocatable :: face_pts(:, :)
142+
143+
allocate (poly(3, 2))
49144

50-
! points on faces
51-
call check(error, point_in_polygon(0.0_DP, 0.5_DP, poly))
52-
call check(error, point_in_polygon(0.5_DP, 0.0_DP, poly))
53-
call check(error, point_in_polygon(1.0_DP, 0.5_DP, poly))
54-
call check(error, point_in_polygon(0.5_DP, 1.0_DP, poly))
145+
allocate (in_pts(3, 2))
146+
in_pts(1, :) = (/0.8_DP, 0.0001_DP/)
147+
in_pts(2, :) = (/0.5_DP, 0.49999_DP/)
148+
in_pts(3, :) = (/0.0001_DP, 0.8_DP/)
149+
150+
allocate (out_pts(2, 2))
151+
out_pts(1, :) = (/0.5_DP, 0.50001_DP/)
152+
out_pts(2, :) = (/-0.5_DP, 34.0_DP/)
153+
154+
allocate (vert_pts(3, 2))
155+
vert_pts(1, :) = (/0.0_DP, 0.0_DP/)
156+
vert_pts(2, :) = (/1.0_DP, 0.0_DP/)
157+
vert_pts(3, :) = (/0.0_DP, 1.0_DP/)
158+
159+
allocate (face_pts(3, 2))
160+
face_pts(1, :) = (/0.0_DP, 0.5_DP/)
161+
face_pts(2, :) = (/0.5_DP, 0.0_DP/)
162+
face_pts(3, :) = (/0.5_DP, 0.5_DP/)
163+
164+
poly(1, :) = (/0.0_DP, 0.0_DP/)
165+
poly(2, :) = (/0.0_DP, 1.0_DP/)
166+
poly(3, :) = (/1.0_DP, 0.0_DP/)
167+
call test_point_in_polygon(error, "clockwise triangle", &
168+
poly, in_pts, out_pts, vert_pts, face_pts)
55169
if (allocated(error)) return
56170

57-
! vertices counter-clockwise
58171
poly(1, :) = (/0.0_DP, 0.0_DP/)
59172
poly(2, :) = (/1.0_DP, 0.0_DP/)
60-
poly(3, :) = (/1.0_DP, 1.0_DP/)
61-
poly(4, :) = (/0.0_DP, 1.0_DP/)
62-
63-
! points inside polygon
64-
call check(error, point_in_polygon(0.99_DP, 0.01_DP, poly))
65-
call check(error, point_in_polygon(0.5_DP, 0.5_DP, poly))
66-
call check(error, point_in_polygon(0.0001_DP, 0.9999_DP, poly))
173+
poly(3, :) = (/0.0_DP, 1.0_DP/)
174+
call test_point_in_polygon(error, "counter-clockwise triangle", &
175+
poly, in_pts, out_pts, vert_pts, face_pts)
67176
if (allocated(error)) return
68177

69-
! points outside polygon
70-
call check(error, (.not. point_in_polygon(0.5_DP, 1.00001_DP, poly)))
71-
call check(error, (.not. point_in_polygon(-0.5_DP, 34.0_DP, poly)))
72-
if (allocated(error)) return
178+
deallocate (poly)
179+
deallocate (in_pts)
180+
deallocate (out_pts)
181+
deallocate (vert_pts)
182+
deallocate (face_pts)
183+
end subroutine test_point_in_polygon_tri
73184

74-
! points on vertices
75-
call check(error, point_in_polygon(0.0_DP, 0.0_DP, poly))
76-
call check(error, point_in_polygon(1.0_DP, 0.0_DP, poly))
77-
call check(error, point_in_polygon(0.0_DP, 1.0_DP, poly))
78-
call check(error, point_in_polygon(1.0_DP, 1.0_DP, poly))
185+
!> @brief Test an irregular polygon
186+
subroutine test_point_in_polygon_irr(error)
187+
type(error_type), allocatable, intent(out) :: error
188+
real(DP), allocatable :: poly(:, :)
189+
real(DP), allocatable :: in_pts(:, :)
190+
real(DP), allocatable :: out_pts(:, :)
191+
real(DP), allocatable :: vert_pts(:, :)
192+
real(DP), allocatable :: face_pts(:, :)
193+
194+
allocate (poly(5, 2))
195+
196+
allocate (in_pts(3, 2))
197+
in_pts(1, :) = (/0.5_DP, 0.1_DP/)
198+
in_pts(2, :) = (/0.5_DP, 0.49_DP/)
199+
in_pts(3, :) = (/1.999_DP, 1.999_DP/)
200+
201+
allocate (out_pts(3, 2))
202+
out_pts(1, :) = (/0.5_DP, -0.1_DP/)
203+
out_pts(2, :) = (/0.5_DP, 0.51_DP/)
204+
out_pts(3, :) = (/-0.5_DP, 34.0_DP/)
205+
206+
allocate (vert_pts(5, 2))
207+
vert_pts(1, :) = (/0.0_DP, 0.0_DP/)
208+
vert_pts(2, :) = (/1.0_DP, 1.0_DP/)
209+
vert_pts(3, :) = (/1.0_DP, 2.0_DP/)
210+
vert_pts(4, :) = (/2.0_DP, 2.0_DP/)
211+
vert_pts(5, :) = (/2.0_DP, 0.0_DP/)
212+
213+
allocate (face_pts(3, 2))
214+
face_pts(1, :) = (/0.5_DP, 0.5_DP/)
215+
face_pts(2, :) = (/2.0_DP, 1.0_DP/)
216+
face_pts(3, :) = (/1.5_DP, 2.0_DP/)
217+
218+
poly(1, :) = (/0.0_DP, 0.0_DP/)
219+
poly(2, :) = (/1.0_DP, 1.0_DP/)
220+
poly(3, :) = (/1.0_DP, 2.0_DP/)
221+
poly(4, :) = (/2.0_DP, 2.0_DP/)
222+
poly(5, :) = (/2.0_DP, 0.0_DP/)
223+
call test_point_in_polygon(error, &
224+
"clockwise irregular polygon", &
225+
poly, in_pts, out_pts, vert_pts, face_pts)
79226
if (allocated(error)) return
80227

81-
! points on faces
82-
call check(error, point_in_polygon(0.0_DP, 0.5_DP, poly))
83-
call check(error, point_in_polygon(0.5_DP, 0.0_DP, poly))
84-
call check(error, point_in_polygon(1.0_DP, 0.5_DP, poly))
85-
call check(error, point_in_polygon(0.5_DP, 1.0_DP, poly))
228+
poly(1, :) = (/0.0_DP, 0.0_DP/)
229+
poly(2, :) = (/2.0_DP, 0.0_DP/)
230+
poly(3, :) = (/2.0_DP, 2.0_DP/)
231+
poly(4, :) = (/1.0_DP, 2.0_DP/)
232+
poly(5, :) = (/1.0_DP, 1.0_DP/)
233+
call test_point_in_polygon(error, &
234+
"counter-clockwise irregular polygon", &
235+
poly, in_pts, out_pts, vert_pts, face_pts)
86236
if (allocated(error)) return
87237

88-
end subroutine test_point_in_polygon_rect
238+
deallocate (poly)
239+
deallocate (in_pts)
240+
deallocate (out_pts)
241+
deallocate (vert_pts)
242+
deallocate (face_pts)
243+
end subroutine test_point_in_polygon_irr
89244

90245
end module TestGeomUtil

‎src/Utilities/GeomUtil.f90

Lines changed: 12 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@ module GeomUtilModule
88
!> @brief Check if a value is between two other values (inclusive).
99
logical function between(x, a, b)
1010
real(DP), intent(in) :: x, a, b
11-
between = (x >= a .and. x <= b .or. x <= a .and. x >= b)
11+
between = ((x >= a .and. x <= b) .or. (x <= a .and. x >= b))
1212
end function between
1313

1414
!> @brief Check if a point is within a polygon.
@@ -23,40 +23,40 @@ logical function point_in_polygon(x, y, poly)
2323

2424
point_in_polygon = .false.
2525
num_verts = size(poly, 1)
26-
xa = poly(1, 1)
27-
ya = poly(1, 2)
26+
xa = poly(num_verts, 1)
27+
ya = poly(num_verts, 2)
2828

29-
do i = 0, num_verts + 1
29+
do i = 0, num_verts - 1
3030
ii = mod(i, num_verts) + 1
3131
xb = poly(ii, 1)
3232
yb = poly(ii, 2)
3333

34-
! boundary cases
3534
if ((x == xa .and. y == ya) .or. (x == xb .and. y == yb)) then
36-
! vertex point
35+
! on vertex
3736
point_in_polygon = .true.
37+
exit
3838
else if (ya == yb .and. y == ya .and. between(x, xa, xb)) then
39-
! horizontal edge
39+
! on horizontal edge
4040
point_in_polygon = .true.
41-
! if within vertical range, cast a ray
41+
exit
4242
else if (between(y, ya, yb)) then
43-
if (y == ya .and. yb >= ya .or. y == yb .and. ya >= yb) then
43+
if ((y == ya .and. yb >= ya) .or. (y == yb .and. ya >= yb)) then
4444
xa = xb
4545
ya = yb
4646
cycle
4747
end if
4848
! cross product
4949
c = (xa - x) * (yb - y) - (xb - x) * (ya - y)
50-
! boundary case
5150
if (c == 0) then
51+
! on edge
5252
point_in_polygon = .true.
53-
! standard intersection
53+
exit
5454
else if ((ya < yb) .eqv. (c > 0)) then
55+
! ray intersection
5556
point_in_polygon = .not. point_in_polygon
5657
end if
5758
end if
5859

59-
if (point_in_polygon) exit
6060
xa = xb
6161
ya = yb
6262
end do

0 commit comments

Comments
 (0)