Calculate the total length of line segments within polygon

Question

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

Any ideas please?

Answer 1

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:

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.

Answer 2

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