Speed Up Sampling of Perimeter of Polygon

I have a polygon whose vertices (points) are defined by some image coordinates that I have manually selected. The polygon is used to create a mask for a region of interest in a microscopy image. Typically I manually select between 10 to 50 points in creating the polygon. I would like to be able to uniformly sample the perimeter of the polygon to get more points along it (with the number of points being some proportion of the total number of points in the perimeter). Below is some code I have been using but it is rather slow and doesn't quite uniformly sample the perimeter. The slow steps seem to be in the creation of the polygonal mask via the Graphics function and the use of FindShortestTour to maintain the ordering of the perimeter points.

samplePolyPerim[pointsListIn_, perimSamplingRatio_] := (

(* create a white polygon on a black background to serve as a mask
to use to designate region of interest on microscopy image of size
600x600 *)

polyMask =
Graphics[{White, Polygon[pointsListIn]},
PlotRange -> {{1, 600}, {1, 600}}, ImageSize -> {600, 600},
AspectRatio -> 1, ImagePadding -> {{0, 0}, {0, 0}},
Background -> Black];

(* find the image coordinates of the perimeter of the polygon.  White
pixels have value of 1. *)

perimPositions =
ImageValuePositions[MorphologicalPerimeter[polyMask], 1];

(* get the proper ordering of the image coordinates to maintain
the shape of the polygon *)

perimPositions$TSPordering = FindShortestTour[perimPositions][[2]]; (* number of points in polygon perimeter to sample is set as a proportion of the number of points found by MorphologicalPerimeter *) totalPerimPointsToSample = Round[perimSamplingRatio*(Length@perimPositions$TSPordering), 1];

(* the delta_sample between index of the point along perimeter *)

samplingDeltaRounded =
Round[Length@perimPositions$TSPordering/totalPerimPointsToSample, 1]; (* taking samples of the FindShortestTour ordering of perimeter points *) perimPositions$TSPordering$Sampled = Table[perimPositions$TSPordering[[i]], {i, 1,
Length@perimPositions$TSPordering, samplingDeltaRounded}]; (* subsetting the perimeter positions from MorphologicalPerimeter based on ordered sampling *) uniformSampledPerimPoints = perimPositions[[perimPositions$TSPordering\$Sampled]];

Return[uniformSampledPerimPoints]

)


Here is an example of passing ~20 manually chosen points to the samplePolyPerim function and it returning ~30 points nearly uniform points along the perimeter (I have used 10% as the proportion of points along the perimeter to sample):

manualPolyPoints = {{355., 396.}, {356., 410.}, {349.,
420.}, {333., 431.}, {323., 432.}, {310., 438.}, {298.,
441.}, {295., 453.}, {293., 464.}, {282., 472.}, {273.,
459.}, {273., 443.}, {272., 433.}, {272., 420.}, {280.,
412.}, {291., 403.}, {308., 398.}, {325., 391.}, {341.,
392.}};

In[13]:= samplePolyPerim[manualPolyPoints, 0.1]

Out[13]= {{281.5, 470.5}, {287.5, 467.5}, {293.5, 462.5}, {294.5,
454.5}, {296.5, 446.5}, {299.5, 439.5}, {307.5, 437.5}, {314.5,
435.5}, {321.5, 431.5}, {330.5, 430.5}, {337.5, 427.5}, {343.5,
422.5}, {349.5, 419.5}, {353.5, 413.5}, {356.5, 407.5}, {355.5,
397.5}, {349.5, 393.5}, {341.5, 391.5}, {332.5, 390.5}, {323.5,
390.5}, {315.5, 393.5}, {308.5, 396.5}, {301.5, 399.5}, {293.5,
401.5}, {287.5, 405.5}, {281.5, 410.5}, {277.5, 414.5}, {272.5,
420.5}, {272.5, 430.5}, {273.5, 440.5}, {273.5, 449.5}, {273.5,
459.5}, {278.5, 465.5}}


My main concern is trying to speed-up the function (AbsoluteTiming is around ~0.20 seconds on my machine) but thoughts on getting better perimeter sampling are welcome too. Also, I am hoping to get some suggestions for finding a faster way to create a mask for my region of interest that I can ultimately use with ComponentMeasurements to get data out of images. Going through Graphics seems slow when I just need a label matrix for the points inside the polygon.

• RandomPoint[RegionBoundary@Polygon[manualPolyPoints], 30]? – Rahul May 8 '16 at 16:37
• @RunnyKine This is definitely faster but is there any way to draw samples that are some equal distance apart? – user13999 May 8 '16 at 16:54
• @Rahul Thanks for this suggestion. My worry with RandomPoint is that the sampling could miss parts of the polygon that are important to the shape of the region. Is there any generic function for a region similar to CirclePoints? – user13999 May 8 '16 at 16:56
• You can use linear interpolation for multisecting a line segment: {1 - Subdivide[n], Subdivide[n]}.{pt1, pt2} – J. M. will be back soon May 8 '16 at 16:58
• @RunnyKine Sorry I don't understand your this comment.You mean this code will doesn't work normally？ – yode May 9 '16 at 6:41

1 Answer

You can use Interpolate to interpolate between the polygon vertices:

First make sure the polygon is closed, by appending the first vertex:

cyclic = Append[manualPolyPoints, First[manualPolyPoints]];


then accumulate the distance from one vertex to the next:

accumulatedDistance =
Rescale@Prepend[Accumulate[Norm /@ Differences[cyclic]], 0.];


then create a linear interpolation from the arc length to x/y coordinates:

{intX, intY} =
Interpolation[Transpose[{accumulatedDistance, #}],
InterpolationOrder -> 1] & /@ Transpose[cyclic];
border = Through[{intX, intY}[#]] &;


Now border is basically arc-length parametrization of your polygon border. You can simply call it with evenly spaced reals to get evenly spaced samples:

samples = border /@ Range[0, 1, 0.01];
ListPlot[samples, AspectRatio -> Automatic,
Prolog -> {Red, Opacity[.3], Line[cyclic]}]


ADD Timing measurement:

uniformSamples[pts_, d_: 0.01] :=
Module[{cyclic, accumulatedDistance, intX, intY, border},
cyclic = Append[pts, First[pts]];
accumulatedDistance =
Rescale@Prepend[Accumulate[Norm /@ Differences[cyclic]], 0.];
{intX, intY} =
Interpolation[Transpose[{accumulatedDistance, #}],
InterpolationOrder -> 1] & /@ Transpose[cyclic];
border = Through[{intX, intY}[#]] &;
border /@ Range[0, 1, d]]

RepeatedTiming[uniformSamples[manualPolyPoints];]


{0.000650, Null}

• Very nice and fast! Thanks. Never used Through before. Do you have any suggestions for how to speedily make a label matrix mask out of the polygon that can be passed to ComponentMeasurements without going through Graphics as I did in the code I posted where I create polyMask? – user13999 May 8 '16 at 18:16
• @user13999: Sorry, no. You can try to ask another question for that. But in general, I'd try to avoid creating a label matrix (unless performance is irrelevant) and calculate the measurement directly from the vertices instead – Niki Estner May 8 '16 at 18:59
• Yes, I will post a new question about this. Not entirely clear on what you mean by "measuring just from the vertices". I understand that can work for spatial stats (area, perimeter, etc) but I also need to extract data on the image pixel intensities. Not sure how to do this without creating a label matrix, which seems to be fairly slow in all the approaches I've tried. – user13999 May 8 '16 at 21:54