Ray-casting algorithm: Difference between revisions

Content added Content deleted
(fortran)
(→‎{{header|Tcl}}: Use the better algorithm from the top of the page to fix the error)
Line 366: Line 366:


=={{header|Tcl}}==
=={{header|Tcl}}==
{{incorrect|Tcl|If you try with {-5 5} {3 0 7 0 10 5 7 10 3 10 0 5 } (exagon), you notice a bug: it says the point is inside. Even if the point is {100 5} indeed! What's wrong?}}
<lang Tcl>package require Tcl 8.5
<lang Tcl>package require Tcl 8.5

proc point_in_polygon {point polygon} {
proc point_in_polygon {point polygon} {
set count 0
set count 0
foreach side [sides $polygon] {
foreach side [sides $polygon] {
if [ray_intersects_line $point $side] {
if {[ray_intersects_line $point $side]} {
incr count
incr count
}
}
Line 379: Line 378:
}
}
proc sides polygon {
proc sides polygon {
foreach {x0 y0} $polygon break
lassign $polygon x0 y0
foreach {x y} [lrange [lappend polygon $x0 $y0] 2 end] {
foreach {x y} [lrange [lappend polygon $x0 $y0] 2 end] {
lappend res [list $x0 $y0 $x $y]
lappend res [list $x0 $y0 $x $y]
Line 388: Line 387:
}
}
proc ray_intersects_line {point line} {
proc ray_intersects_line {point line} {
foreach {x y} $point break
lassign $point Px Py
foreach {xa xb ya yb} $line break
lassign $line Ax Ay Bx By
# Reverse line direction if necessary
set xout [expr {max($xa,$xb) + 1}]
lines_cross [list $x $y $xout $y] $line
if {$By < $Ay} {
lassign $line Bx By Ax Ay
}
}
proc lines_cross {line1 line2} {
# Add epsilon to
foreach {xa ya xb yb} $line1 break
foreach {xc yc xd yd} $line2 break
if {$Py == $Ay || $Py == $By} {
set Py [expr {$Py + abs($Py)/1e6}]

}
set det [expr {($xb - $xa) * ($yd - $yc) - ($xd - $xc) * ($yb - $ya)}]
# Bounding box checks
if {$det == 0} {return 0} ;#-- parallel or collinear
if {$Py < $Ay || $Py > $By || $Px > max($Ax,$Bx)} {

return 0
set nump [expr {($xd - $xc) * ($yb + $ya) - ($xb + $xa) * ($yd - $yc) + 2. * ($xc * $yd - $xd * $yc)}]
} elseif {$Px < min($Ax,$Bx)} {
if {abs($nump) > abs($det)} {return 0} ;#-- cross outside line1
return 1

}
set numq [expr {($xd + $xc) * ($yb - $ya) - ($xb - $xa) * ($yd + $yc) + 2. * ($xb * $ya - $xa * $yb)}]
# Compare dot products to compare (cosines of) angles
if {abs($numq) > abs($det)} {return 0} ;#-- cross outside line2
set mRed [expr {$Ax != $Bx ? ($By-$Ay)/($Bx-$Ax) : Inf}]

set mBlu [expr {$Ax != $Px ? ($Py-$Ay)/($Px-$Ax) : Inf}]
return 1
return [expr {$mBlu >= $mRed}]
}
}

foreach {point poly} {
foreach {point poly} {
{0 0} {-1 -1 -1 1 1 1 1 -1}
{0 0} {-1 -1 -1 1 1 1 1 -1}
Line 420: Line 420:
{10 10} {0 0 2.5 2.5 0 10 2.5 7.5 7.5 7.5 10 10 10 0 7.5 0.1}
{10 10} {0 0 2.5 2.5 0 10 2.5 7.5 7.5 7.5 10 10 10 0 7.5 0.1}
{2.5 2.5} {0 0 2.5 2.5 0 10 2.5 7.5 7.5 7.5 10 10 10 0 7.5 0.1}
{2.5 2.5} {0 0 2.5 2.5 0 10 2.5 7.5 7.5 7.5 10 10 10 0 7.5 0.1}
{-5 5} {3 0 7 0 10 5 7 10 3 10 0 5 }
} {
} {
puts "$point in $poly = [point_in_polygon $point $poly]"
puts "$point in $poly = [point_in_polygon $point $poly]"