3
$\begingroup$

So we have a polygon with N vertices located on grid. All vertices are located at the intersection of cells (so their coordinates are integers).

The objective is to calculate the total length of line segments within given polygon. There is the simple way to do this when given polygon is right triangle.

We can just find the length of lines inside rectangle which contains our triangle

xSegments = x(y - 1)
ySegments = y(x - 1)
rectangleSegments = xSegments + ySegments = 2 x y - (x + y)
triangleSegments = totalSegments / 2 = x y - (x + y)/2

But what should we do with other types of polygon like this:

Graphics[
    {
    EdgeForm[Black], FaceForm[None],
    Polygon[{{0,0}, {4,-1}, {4,1}, {2,1}, {0,3}}],
    },
    GridLines->Automatic,
    PlotRange->{{-1,5}, {-2,4}}
]

enter image description here

Any ideas please?

$\endgroup$
  • $\begingroup$ If I understood you correctly: you want the perimeter of your lattice polygons? $\endgroup$ – J. M. will be back soon May 2 '15 at 10:55
  • $\begingroup$ I don't actually know what is perimeter of lattice polygons. Need to find count of all grey lines inside the polygon (see pictures). In my task they are called line segments. Look at the picture. I am talking about green lines: s4.postimg.org/kgpm37z8t/image.png $\endgroup$ – Vladyslav May 2 '15 at 11:00
  • $\begingroup$ If this question is really about the software Mathematica, then you should correct the syntax of your code. xy is not the same as y*x, which is equivalent to y x. $\endgroup$ – Karsten 7. May 2 '15 at 11:06
5
$\begingroup$

supposing that your polynom is given with 2 dimensional points in counter-clockwise ordering. For example:

points = {{4, 0}, {1, 2}, {0, 4}, {-3, 0}, {0, -1}};

Lets define a "test area" in which the polynom is embeded:

{minX, maxX} = {Min@#, Max@#} &@points〚All, 1〛;
{minY, maxY} = {Min@#, Max@#} &@points〚All, 2〛;

The edges of the polynom are lines which are defined with 2 corner points:

lines = Partition[points, 2, 1, 1];

The following will calculate the intersection points with the horizontal lines, based on basical analysis (y = M x + B)

intersectionPointsX = {};
For[y = minY + 1, y < maxY, y++,

crossingLines = Select[lines, #[[1, 2]] =!= #[[2, 2]] &];

For[j = 1, j <= Length[crossingLines], j++,
 {ux, uy} = crossingLines〚j, 1〛;
 {vx, vy} = crossingLines〚j, 2〛;
 If[Min[uy, vy] <= y <= Max[uy, vy],
   If[ux == vx,
     AppendTo[intersectionPointsX, {ux, y}],

     M = -((-uy + vy)/(ux - vx)); B = -((uy vx - ux vy)/(ux - vx));
     AppendTo[intersectionPointsX, {(y - B)/M, y}]
   ]
 ]
]
]

Lets do the same with the vertical lines:

intersectionPointsY = {};
For[x = minX + 1, x <  maxX, x++,

 crossingLines = Select[lines, #[[1, 1]] =!= #[[2, 1]] &];

 For[j = 1, j <= Length[crossingLines], j++,
  {ux, uy} = crossingLines〚j, 1〛;
  {vx, vy} = crossingLines〚j, 2〛;
  If[Min[ux, vx] <= x <= Max[ux, vx],
    If[ux == vx,
    AppendTo[intersectionPointsY, {x, (uy+vy)/2}]
    , 
    M = -((-uy + vy)/(ux - vx)); B = -((uy vx - ux vy)/(ux - vx));
    AppendTo[intersectionPointsY, {x, M x + B}]
    ]
  ]
 ]
]

Now we Gather these intersection points with respect to their intersecting line. We also delete possible dublicates.

intersectionPointsX = 
DeleteDuplicates /@ GatherBy[intersectionPointsX, Last];
intersectionPointsY = 
DeleteDuplicates /@ GatherBy[intersectionPointsY, First];

Edit: The following transformations will include more special cases, if there are more than 2 intersections for one line:

intersectionPointsX = 
 If[Length[#] > 2, 
   Sort[Select[Partition[#, 2, 1], ! MemberQ[lines, #] &], 
   Abs[#1[[1, 1]] - #1[[2, 1]]] < Abs[#2[[1, 1]] - #2[[2, 1]]] &][[
  1]], #] & /@ intersectionPointsX

intersectionPointsY = 
 If[Length[#] > 2, 
   Sort[Select[Partition[#, 2, 1], ! MemberQ[lines, #] &], 
   Abs[#1[[1, 1]] - #1[[2, 1]]] < 
     Abs[#2[[1, 1]] - #2[[2, 1]]] &][[1]], #] & /@ 
 intersectionPointsY;

Lets calculate the total segmentslenth in x- and y direction:

ΔX = 
Total[Abs[Subtract @@@ intersectionPointsX〚All, All, 1〛]] // N;
ΔY = 
Total[Abs[Subtract @@@ intersectionPointsY〚All, All, 2〛]] // N;

so the total lenth of line segments is:

totalSegmentsLenth = ΔX + ΔY;

In this example it is:

31

Visualization

ListPlot[{points, Flatten[intersectionPointsX, 1], 
          Flatten[intersectionPointsY, 1]}, 

         PlotStyle -> {Directive[Red, PointSize[0.02], Opacity[.5]], 
           Directive[Blue, PointSize[0.02], Opacity[.5]], 
           Directive[Green, PointSize[0.02]]}, 

         Prolog -> {Line[points~Join~{points〚1〛}], Opacity[.1], 
           Polygon[points]},

         Epilog -> {Text[#1, Offset[{10, 10}, #2]] & @@@ 
           Transpose[{FromCharacterCode /@ 
           Range[65, 65 + Length[points] - 1], points}],
           Thickness[0.005],
           Blue, Line /@ intersectionPointsX, Green, 
           Line /@ intersectionPointsY},

        PlotRangePadding -> Scaled[.2],
        Frame -> True,
        GridLinesStyle -> GrayLevel[.8],
        GridLines -> {Table[x, {x, minX, maxX}],
          Table[y, {y, minY,maxY}]}, 
        Axes -> False, 
        PlotLabel -> Column[{
          Style["ΔX = "<>ToString@ΔX, Blue],
          Style["ΔY = "<>ToString@ΔY,Darker@Green], 
          "ΔTotal = "<>ToString[totalSegmentsLenth]}]
]

Out:

enter image description here

Another example calculated with the same code:

points = {{1, 2}, {-1, 2}, {-1, -2}, {1, -2}};

enter image description here

One can see that horizontal and vertical edges are not counted as they where inside the polygon. You can include these edges by using For[y = minY, y <= maxY, ...] instead of For[y = minY + 1, y < maxY, ...] and the same with x.

Here is an example that produced errors (division by zero, see comments) with the previous version of this code:

points = {{0, 0}, {0, 2}, {-2, 2}, {-4, 4}, {-4, 1}, {-5, -1}, {-3, 
0}};

enter image description here

non convex polygons may not work properly:

points = {{0, 0}, {0, 4}, {-2, 2}, {-4, 4}, {-4, 1}, {-4, 0}, {-3, 0}};

enter image description here

You could cut a non convex polygon into pieces of convex polygons.

$\endgroup$
  • $\begingroup$ I am going to use only a part of this answer, but it's better to correct it for others. Your code which begins with 'Prolog' and 'Frame' isn't in code tag so it's hard to read. $\endgroup$ – Vladyslav May 2 '15 at 17:13
  • $\begingroup$ @Seprum I needed some time for formatting. The code is relatively huge. But thanks for the hint. $\endgroup$ – sacratus May 2 '15 at 17:28
  • $\begingroup$ I will try to transfer the algorithm's code into Python later, so if there will be some questions about it or about syntax, I will ask you, okay? $\endgroup$ – Vladyslav May 2 '15 at 17:33
  • $\begingroup$ @Seprum I have no experience with python, but you can ask if you have problems with understanding the mathematica syntax. $\endgroup$ – sacratus May 2 '15 at 17:41
  • $\begingroup$ In my polygon (at the picture) there is devision by zero error attempting to initialize M = -((-uy + vy) / (ux - vx)), because the line is vertical. What am I supposed to do in that case? And why are u checking if ux is equal to vx? The division can't be done in that case I think. $\endgroup$ – Vladyslav May 3 '15 at 8:32
1
$\begingroup$

I think the simplest approach is to use Region functionality. Create possible infinite lines that might intersect with your region, and then find the measure of the intersection of your region with the lines. Here is a function that does this:

XYSegmentLength[region_] := With[{bounds = RegionBounds[region]},
    RegionMeasure @ RegionIntersection[
        region,
        RegionUnion @@ Join[
            Table[
                InfiniteLine[{{x, 0}, {x, 1}}],
                {x, Ceiling[bounds[[1, 1]]], Floor[bounds[[1,-1]]]}
            ],
            Table[
                InfiniteLine[{{0, x}, {1, x}}],
                {x, Ceiling[bounds[[2, 1]]], Floor[bounds[[2, -1]]]}
            ]
        ]
    ]
]

For your example:

XYSegmentLength[Polygon[{{0,0}, {4,-1}, {4,1}, {2,1}, {0,3}}]]

39/2

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.