# Calculate the total length of line segments within polygon

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}}
] • If I understood you correctly: you want the perimeter of your lattice polygons? May 2, 2015 at 10:55
• 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 May 2, 2015 at 11:00
• 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. May 2, 2015 at 11:06

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]]] &][], #] & /@
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},

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: Another example calculated with the same code:

points = {{1, 2}, {-1, 2}, {-1, -2}, {1, -2}}; 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}}; non convex polygons may not work properly:

points = {{0, 0}, {0, 4}, {-2, 2}, {-4, 4}, {-4, 1}, {-4, 0}, {-3, 0}}; You could cut a non convex polygon into pieces of convex polygons.

• 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. May 2, 2015 at 17:13
• @Seprum I needed some time for formatting. The code is relatively huge. But thanks for the hint. May 2, 2015 at 17:28
• 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? May 2, 2015 at 17:33
• @Seprum I have no experience with python, but you can ask if you have problems with understanding the mathematica syntax. May 2, 2015 at 17:41
• 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. May 3, 2015 at 8:32

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]]]}
]
]
]
]


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