Anonymous user
Ray-casting algorithm: Difference between revisions
fortran
m (inverted the logic when modified A below B and viceversa; fixed) |
(fortran) |
||
Line 198:
The test's output reveals the meaning of <cite>coherent results</cite>: a point on the leftmost vertical side of the square (coordinate 0,5) is considered outside; while a point on the rightmost vertical side of the square (coordinate 10,5) is considered inside, but on the top-right vertex (coordinate 10,10), it is considered outside again.
=={{header|Fortran}}==
{{works with|Fortran|95 and later}}
The following code uses the <tt>Points_Module</tt> defined [[Closest pair problem#Fortran|here]].
This module ''defines'' "polygons".
<lang fortran>module Polygons
use Points_Module
implicit none
type polygon
type(point), dimension(:), allocatable :: points
integer, dimension(:), allocatable :: vertices
end type polygon
contains
function create_polygon(pts, vt)
type(polygon) :: create_polygon
type(point), dimension(:), intent(in) :: pts
integer, dimension(:), intent(in) :: vt
integer :: np, nv
np = size(pts,1)
nv = size(vt,1)
allocate(create_polygon%points(np), create_polygon%vertices(nv))
create_polygon%points = pts
create_polygon%vertices = vt
end function create_polygon
subroutine free_polygon(pol)
type(polygon), intent(inout) :: pol
deallocate(pol%points, pol%vertices)
end subroutine free_polygon
end module Polygons</lang>
The ray casting algorithm module:
<lang fortran>module Ray_Casting_Algo
use Polygons
implicit none
real, parameter, private :: eps = 0.00001
private :: ray_intersects_seg
contains
function ray_intersects_seg(p0, a0, b0) result(intersect)
type(point), intent(in) :: p0, a0, b0
logical :: intersect
type(point) :: a, b, p
real :: m_red, m_blue
p = p0
! let variable "a" be the point with smallest y coordinate
if ( a0%y > b0%y ) then
b = a0
a = b0
else
a = a0
b = b0
end if
if ( (p%y == a%y) .or. (p%y == b%y) ) p%y = p%y + eps
intersect = .false.
if ( (p%y > b%y) .or. (p%y < a%y) ) return
if ( p%x > max(a%x, b%x) ) return
if ( p%x < min(a%x, b%x) ) then
intersect = .true.
else
if ( abs(a%x - b%x) > tiny(a%x) ) then
m_red = (b%y - a%y) / (b%x - a%x)
else
m_red = huge(m_red)
end if
if ( abs(a%x - p%x) > tiny(a%x) ) then
m_blue = (p%y - a%y) / (p%x - a%x)
else
m_blue = huge(m_blue)
end if
if ( m_blue >= m_red ) then
intersect = .true.
else
intersect = .false.
end if
end if
end function ray_intersects_seg
function point_is_inside(p, pol) result(inside)
logical :: inside
type(point), intent(in) :: p
type(polygon), intent(in) :: pol
integer :: i, cnt, pa, pb
cnt = 0
do i = lbound(pol%vertices,1), ubound(pol%vertices,1), 2
pa = pol%vertices(i)
pb = pol%vertices(i+1)
if ( ray_intersects_seg(p, pol%points(pa), pol%points(pb)) ) cnt = cnt + 1
end do
inside = .true.
if ( mod(cnt, 2) == 0 ) then
inside = .false.
end if
end function point_is_inside
end module Ray_Casting_Algo</lang>
'''Testing'''
<lang fortran>program Pointpoly
use Points_Module
use Ray_Casting_Algo
implicit none
character(len=16), dimension(4) :: names
type(polygon), dimension(4) :: polys
type(point), dimension(14) :: pts
type(point), dimension(7) :: p
integer :: i, j
pts = (/ point(0,0), point(10,0), point(10,10), point(0,10), &
point(2.5,2.5), point(7.5,2.5), point(7.5,7.5), point(2.5,7.5), &
point(0,5), point(10,5), &
point(3,0), point(7,0), point(7,10), point(3,10) /)
polys(1) = create_polygon(pts, (/ 1,2, 2,3, 3,4, 4,1 /) )
polys(2) = create_polygon(pts, (/ 1,2, 2,3, 3,4, 4,1, 5,6, 6,7, 7,8, 8,5 /) )
polys(3) = create_polygon(pts, (/ 1,5, 5,4, 4,8, 8,7, 7,3, 3,2, 2,5 /) )
polys(4) = create_polygon(pts, (/ 11,12, 12,10, 10,13, 13,14, 14,9, 9,11 /) )
names = (/ "square", "square hole", "strange", "exagon" /)
p = (/ point(5,5), point(5, 8), point(-10, 5), point(0,5), point(10,5), &
point(8,5), point(10,10) /)
do j = 1, size(p)
do i = 1, size(polys)
write(*, "('point (',F8.2,',',F8.2,') is inside ',A,'? ', L)") &
p(j)%x, p(j)%y, names(i), point_is_inside(p(j), polys(i))
end do
print *, ""
end do
do i = 1, size(polys)
call free_polygon(polys(i))
end do
end program Pointpoly</lang>
=={{header|Tcl}}==
| |||