MathGroup Archive 2001

[Date Index] [Thread Index] [Author Index]

Search the Archive

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>
  • Prev by Date: Re: 2 questions about bitmap images
  • Next by Date: Re: list of bits to string
  • Previous by thread: Re: Interior of a polygon
  • Next by thread: Plotting