Ray-casting algorithm: Difference between revisions
Content added Content deleted
m (inverted the logic when modified A below B and viceversa; fixed) |
(fortran) |
||
Line 198: | 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. |
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}}== |
=={{header|Tcl}}== |