Re: Interior of a polygon
- To: mathgroup at smc.vnet.net
- Subject: [mg28686] Re: Interior of a polygon
- From: "Allan Hayes" <hay at haystack.demon.co.uk>
- Date: Wed, 9 May 2001 04:20:07 -0400 (EDT)
- References: <9d85uu$862@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Marius, Andrzej,
Here is a variant of my earlier technique; it generate the filling points
directly from the given simple polygon without triangulation and without
needing to digitize the polygon perimeter.
The option IncludeBorder lets the user choose to include or exclude
integer points on the border.
Options[fillPolygon]={IncludeBorder ->True};
fillPolygon[v_, (opts___)?OptionQ] :=
Module[{ib, segments, exposedSegs, lineToPoints, toPoints, vp,
vpc, slices, pairs, fill},
ib = IncludeBorder /. Flatten[{opts}] /. Options[fillPolygon];
segments = Split[v, #1[[1]] === #2[[1]] & ];
If[ib, segments = Replace[segments, {{x_, a_}, ___, {x_, b_}} :>
Thread[{x, Range[a, b, Sign[b - a]]}], {-3}];
];
exposedSegs =
Cases[Partition[segments, 3, 1, {2, 2}],
{{{a_, _}, ___}, p:{{b_, _}, ___}, {{c_, _}, ___}} /;
(a - b)*(c - b) > 0 -> p
];
lineToPoints[{p:{___, {a_, b_}}, {{c_, d_}, ___}}] /; a =!= c :=
List /@
Prepend[
Transpose[
{#1,
Interpolation[{{a, b}, {c, d}},
InterpolationOrder -> 1] /@ #1
} & [
Take[Range[a, c, Sign[c - a]], {2, -2}]
]
],
p
];
toPoints[u_] :=
Level[lineToPoints /@ Partition[u, 2, 1, {1, 1}], {-3}];
vp = toPoints[segments];
vpc = Complement[vp, exposedSegs];
vpc = Sort[vpc, #1[[1,1]] <= #2[[1,1]] & ];
vpc = Map[Sort, vpc, {-3}];
slices = Split[vpc, #1[[1,1]] == #2[[1,1]] & ];
pairs = (Partition[#1, 2] & ) /@ slices;
fill =
Replace[pairs,
{{___, {x_, a_}}, {{x_, b_}, ___}} :>
Thread[{x,
Range[a /. {u_Integer /; !ib -> u + 1, a -> Ceiling[a]},
b /. {u_Integer /; !ib -> u - 1, b -> Floor[b]}
]}
],
{-4}
];
If[ib, Union, Complement][
DeleteCases[Level[fill, {-2}], {}], Level[segments, {-2}]
]
]
EXAMPLE
v={{10,10},{13,15},{10,30},{15,20},{30,30},{30,28},{30,25},{35,24},{35,
20},{30,17},{30,14},{20,17},{25,5},{15,12}};
shv=Show[Graphics[Line[Append[v,v[[1]]]]],AspectRatio\[Rule]Automatic];
Show[Graphics[{PointSize[.02],Hue[0],Point/@pi}],shv,
AspectRatio\[Rule]Automatic];
Developer`ClearCache[];
pi=fillPolygon[v,IncludeBorder\[Rule]False];//Timing
{0.17 Second,Null}
--
Allan
---------------------
Allan Hayes
Mathematica Training and Consulting
Leicester UK
www.haystack.demon.co.uk
hay at haystack.demon.co.uk
Voice: +44 (0)116 271 4198
Fax: +44 (0)870 164 0565
"Andrzej Kozlowski" <andrzej at tuins.ac.jp> wrote in message
news:9d85uu$862 at smc.vnet.net...
> I have found a little time to implement the solution to Mariusz
Jankowski's
> problem which I sketched in an earlier posting in reply to his original
> message. It only works for convex polygons but can be easily extended to
all
> by means of Martin Kraus' Polygonal trianulation package. The main be bugs
> because I have tested only for a single polygon (given below). The only
> improvement on the original idea which I described is to use the
> GenericCylindricalAlgebraicDecomposition instead of
> CylindricalAlgebraicDecomposition, which ought to make the function much
> faster.
>
>
>
> Here is my function:
>
> <<<< Experimental`
>
> PointsInsidePolygon[poly_List] :=
> Module[{ineq,
> rule1 = {Less[a_, x_, b_] :> {x, Ceiling[a], Floor[b]},
> Inequality[a_, Less, x_, Less, b_] :> {x, Ceiling[a],
Floor[b]}},
> rule2 = {{u_, w_} :> Table[{x, y}, u, w]}},
> ineq[{{a_, b_}, {a_, d_}, {e_, f_}}] := If[e > a, x > a, x < a];
> ineq[{{a_, b_}, {c_, d_}, {e_, f_}}] :=
> If[f > ((d - b)/(c - a))(e - a) + b, y > ((d - b)/(c - a))(x - a) +
b,
> y < ((d - b)/(c - a))(x - a) + b];
> Complement[
> Join @@ Join @@ (((List @@@
> First[GenericCylindricalAlgebraicDecomposition[
> Map[ineq, Partition[poly, 3, 1, {1, 1}]], {x,
> y}]]) /. rule1) /. rule2), poly]]
>
> Note that the argument is a list of vertices in which the first and last
are
> not the same (in other words you do not repeat the starting vetex).
>
> Here is an example:
>
>
>
> polyg = {{1, 3}, {2, 5}, {4, 4}, {5, 2}};
>
> In[55]:=
> pts=PointsInsidePolygon[polyg]
>
> Out[55]=
> {{2,3},{2,4},{3,3},{3,4},{4,3}}
>
> In[56]:=
> s1=Graphics[{RGBColor[1,0,0],Polygon[polyg]}];
>
> In[57]:=
> s2=Graphics[{PointSize[0.02],Point[#]}&/@pts];
>
> In[58]:=
> Show[s1,s2]
>
>
> Andrzej Kozlowski
> Toyama International University
> JAPAN
>
> http://platon.c.u-tokyo.ac.jp/andrzej/
> http://sigma.tuins.ac.jp/~andrzej/
>
>
- Follow-Ups:
- Plotting
- From: Yahsi <yahsi@marun.edu.tr>
- Plotting