PointWithinShape
A set of functions for testing if a point lies within an area
- PointWithinShape(shape, x, y)
- returns true if (x,y) lays within the bounds of the given shape. Shape defined by an array of 0 to n {x=n, y=y} values.
- It generally defers to the other other functions in the file to do the heavy lifting (CrossingsMultipyTest(...), PointWithinLine(...), BoundingBox(...) and colinear(...))
- BoundingBox(box, x, y)
- returns true if the given point is inside the box defined by two points (in the same format taken by PointWithinShape(...)).
- colinear(line, x, y, e)
- returns true if the given point lies on the infinite line defined by two points (again, in the same format taken by PointWithinShape(...))
- e is optional, if given it controls how close (x,y) must be to the line to be considered colinear. it defaults to 0.1
- PointWithinLine(line, x, y, e)
- returns true if the given point lies on the finite line defined by two points (again, in the same format taken by PointWithinShape(...))
- e is optional, if given it controls how close (x,y) must be to the line to be considered within the line. it defaults to 0.66
- CrossingsMultiplyTest(polygon, x, y)
- returns true if the given point lies within the area of the polygon defined by three or points (again, in the same format taken by PointWithinShape(...))
- (it is based on code from Point in Polygon Strategies)
- GetIntersect( points )
- returns the x,y coordinates of the intersect, same format taken by PointWithinShape(...)
- (Added by Bambo - You will need the checkIntersect from the "General math" snippet for it to work)
function PointWithinShape(shape, tx, ty)
if #shape == 0 then
return false
elseif #shape == 1 then
return shape[1].x == tx and shape[1].y == ty
elseif #shape == 2 then
return PointWithinLine(shape, tx, ty)
else
return CrossingsMultiplyTest(shape, tx, ty)
end
end
function BoundingBox(box, tx, ty)
return (box[2].x >= tx and box[2].y >= ty)
and (box[1].x <= tx and box[1].y <= ty)
or (box[1].x >= tx and box[2].y >= ty)
and (box[2].x <= tx and box[1].y <= ty)
end
function colinear(line, x, y, e)
e = e or 0.1
m = (line[2].y - line[1].y) / (line[2].x - line[1].x)
local function f(x) return line[1].y + m*(x - line[1].x) end
return math.abs(y - f(x)) <= e
end
function PointWithinLine(line, tx, ty, e)
e = e or 0.66
if BoundingBox(line, tx, ty) then
return colinear(line, tx, ty, e)
else
return false
end
end
-------------------------------------------------------------------------
-- The following function is based off code from
-- [ http://erich.realtimerendering.com/ptinpoly/ ]
--
--[[
======= Crossings Multiply algorithm ===================================
* This version is usually somewhat faster than the original published in
* Graphics Gems IV; by turning the division for testing the X axis crossing
* into a tricky multiplication test this part of the test became faster,
* which had the additional effect of making the test for "both to left or
* both to right" a bit slower for triangles than simply computing the
* intersection each time. The main increase is in triangle testing speed,
* which was about 15% faster; all other polygon complexities were pretty much
* the same as before. On machines where division is very expensive (not the
* case on the HP 9000 series on which I tested) this test should be much
* faster overall than the old code. Your mileage may (in fact, will) vary,
* depending on the machine and the test data, but in general I believe this
* code is both shorter and faster. This test was inspired by unpublished
* Graphics Gems submitted by Joseph Samosky and Mark Haigh-Hutchinson.
* Related work by Samosky is in:
*
* Samosky, Joseph, "SectionView: A system for interactively specifying and
* visualizing sections through three-dimensional medical image data",
* M.S. Thesis, Department of Electrical Engineering and Computer Science,
* Massachusetts Institute of Technology, 1993.
*
--]]
--[[ Shoot a test ray along +X axis. The strategy is to compare vertex Y values
* to the testing point's Y and quickly discard edges which are entirely to one
* side of the test ray. Note that CONVEX and WINDING code can be added as
* for the CrossingsTest() code; it is left out here for clarity.
*
* Input 2D polygon _pgon_ with _numverts_ number of vertices and test point
* _point_, returns 1 if inside, 0 if outside.
--]]
function CrossingsMultiplyTest(pgon, tx, ty)
local i, yflag0, yflag1, inside_flag
local vtx0, vtx1
local numverts = #pgon
vtx0 = pgon[numverts]
vtx1 = pgon[1]
-- get test bit for above/below X axis
yflag0 = ( vtx0.y >= ty )
inside_flag = false
for i=2,numverts+1 do
yflag1 = ( vtx1.y >= ty )
--[[ Check if endpoints straddle (are on opposite sides) of X axis
* (i.e. the Y's differ); if so, +X ray could intersect this edge.
* The old test also checked whether the endpoints are both to the
* right or to the left of the test point. However, given the faster
* intersection point computation used below, this test was found to
* be a break-even proposition for most polygons and a loser for
* triangles (where 50% or more of the edges which survive this test
* will cross quadrants and so have to have the X intersection computed
* anyway). I credit Joseph Samosky with inspiring me to try dropping
* the "both left or both right" part of my code.
--]]
if ( yflag0 ~= yflag1 ) then
--[[ Check intersection of pgon segment with +X ray.
* Note if >= point's X; if so, the ray hits it.
* The division operation is avoided for the ">=" test by checking
* the sign of the first vertex wrto the test point; idea inspired
* by Joseph Samosky's and Mark Haigh-Hutchinson's different
* polygon inclusion tests.
--]]
if ( ((vtx1.y - ty) * (vtx0.x - vtx1.x) >= (vtx1.x - tx) * (vtx0.y - vtx1.y)) == yflag1 ) then
inside_flag = not inside_flag
end
end
-- Move to the next pair of vertices, retaining info as possible.
yflag0 = yflag1
vtx0 = vtx1
vtx1 = pgon[i]
end
return inside_flag
end
function GetIntersect( points )
local g1 = points[1].x
local h1 = points[1].y
local g2 = points[2].x
local h2 = points[2].y
local i1 = points[3].x
local j1 = points[3].y
local i2 = points[4].x
local j2 = points[4].y
local xk = 0
local yk = 0
if checkIntersect({x=g1, y=h1}, {x=g2, y=h2}, {x=i1, y=j1}, {x=i2, y=j2}) then
local a = h2-h1
local b = (g2-g1)
local v = ((h2-h1)*g1) - ((g2-g1)*h1)
local d = i2-i1
local c = (j2-j1)
local w = ((j2-j1)*i1) - ((i2-i1)*j1)
xk = (1/((a*d)-(b*c))) * ((d*v)-(b*w))
yk = (-1/((a*d)-(b*c))) * ((a*w)-(c*v))
else
xk,yk = 0,0
end
return xk, yk
end