# Finding the corners of noisy polygons

I have some polygons that look for example like this:

If I zoom in very close on one side, you can see the noise.

The data is a list of x coordinates and a corresponding list of y coordinates.

I want an algorithm that will find a much smaller, simpler, less noisy list of coordinates.

I figured that the sides of the polygon are a sequential list of linear equations.

from sklearn.linear_model import Lasso
import numpy as np
import matplotlib.pyplot as plt

xs_name = "xs.txt"
ys_name = "ys.txt"

reg = Lasso(alpha=0.1)
reg.fit(xs, ys)


For reference the xs and ys look like this:

xs

ys

However I only get one coefficient

reg.coef_
Out  [31]: array([0.82647029])


I expected to get a list of coefficients, for each line identified.

I feel like I have conceptually missed something. I'm not even sure Lasso is the right tool for this job.

Does anyone know how I can correctly use Lasso, or alternatively point me to the correct tool for the job.

EDIT

x coords

y coords

I also thought it was worth mentioning that group lasso also got a mention in the survey paper for the ruptures library

EDIT

To be clear, the zoomed in area is circled in red

EDIT

I've had a little success trying an autoregressive model in R

xs <- read.table('xs.txt', sep="\n")

xs <- as.numeric(as.character(unlist(xs)))
ys <- as.numeric(as.character(unlist(ys)))

fastcp_xs <- fastcpd::fastcpd.ar(xs, 3, r.progress = FALSE)
summary(fastcp_xs)
plot(fastcp_xs)


However it seems like the success of this approach may have been mostly luck in this case, as trying this on more data revealed bad results.

Trying the same method on the ys:

fastcp_ys <- fastcpd::fastcpd.ar(ys, 3, r.progress = FALSE)
summary(fastcp_ys)
plot(fastcp_ys)


The autoregressive model was unable to detect the edges for the ys.

The other routines in the fastcpd library seemed to give similarly bad results in my case.

I'm currently thinking my best bet is some form of lasso algorithm. Since the concept of lasso is to fit a sequence of straight lines.

This may turn into a linear programming problem. Maybe I will need to resort to using something like pyomo?

• Interesting question. I don't have anything to contribute, but you might get better answers if you could post the original noisy polygon data somewhere and link to it here. Commented May 16 at 8:52
• not too sure that you can do it with lasso, but most likely you can rely on "old" CV techniques, in particular you can smooth the two coordinates with a gaussian filter, and then check the derivative of neighbors points, that way you can detect some notion of corners cs.cornell.edu/courses/cs664/2003fa/handouts/… Commented May 16 at 11:09
• I noticed lasso was one of the methods used here github.com/doccstat/fastcpd?tab=readme-ov-file#examples
– sav
Commented May 16 at 23:47
• Intuitively, I feel like this is much more difficult for extremely non-convex polygons. Maybe check out some of the literature on statistical shape analysis? Commented May 17 at 0:38
• There are a great many ways to accomplish this. They depend partly on what the polygons represent, on what auxiliary data you might have, on the nature of the noise, and on your accuracy needs. As one simple example of a standard algorithm look up Douglas-Poicker (or Douglas-Peucker as it was originally called). There is a host of methods available through image processing: convert your polygon to a high-res image, use morphological operators to simplify it, and convert back. Etc., etc.
– whuber
Commented May 20 at 2:40

The following (using Mathematica) does not do what one would do "by eye" but that's because the data points don't fall perfectly on a desired number of line segments.

This uses Mathematica's ConcaveHullMesh function with a tuning parameter $$\alpha$$. These "concave hull" shapes are call "alpha shapes". I don't know if the algorithm suggested by @whuber is what ConcaveHullMesh uses.

The x and y files are read into variables x and y (which I don't show). But starting with those variables I "shift" the data by subtracted a number close to the minimum of those variables.

data = Transpose[{x - 225000, y - 6615000}];
α =.;
p = ConstantArray[0, 12];
Do[poly = MeshCoordinates[ConcaveHullMesh[data, i]];
poly = Join[poly, {poly[[1]]}];
p[[i]] = ListPlot[{data, poly, poly}, Joined -> {True, True, False},
PlotStyle -> {LightGray, Blue, {PointSize[0.01], Red}},
ImageSize -> Large, Frame -> True,
PlotLabel -> Style["α = " <> ToString[i], 24, Bold],
PlotLegends -> {"Data points connected", "Convex hull boundary",
"Convex hull points"}],
{i, {1, 2, 3, 4, 5, 6, 7}}]
p[[8]] = p[[6]];
p[[9]] = p[[5]];
p[[10]] = p[[4]];
p[[11]] = p[[3]];
p[[12]] = p[[2]];
Export["poly.gif", p, "DisplayDurations" -> 1]


One can see that with the tuning parameter of 6 there is a "visually adequate" approximation with many fewer points than the original dataset. (This only takes about half a second to complete the polygons and about 10 seconds to create the animated gif.)

• I like that you started this by transposing the coordinates. I might try something like that as it can sometimes result in better numerical stability. The data files I'm using were actually created from an alpha shape. The Douglas-Peucker algorithm suggested by whuber is definitely a different algorithm to the alpha shape algorithm. So, now I have used both of these algorithms together. If I crank up the alpha value too much to get rid of the noise, I lose too much or the 'real' shape. As shown in your animation. The Douglas-Peucker algorithm seems to do well here though.
– sav
Commented May 22 at 6:44
• The answer by C. E. (mathematica.stackexchange.com/questions/167108/…) uses the Douglas-Peucker algorithm. For your example it takes about a minute to complete. However, that can be considerably reduced in time by first using ConcaveHullMesh followed by the Douglas-Peucker algorithm.
– JimB
Commented May 22 at 14:57

Here's a shot at it using the Douglas-Peuker algorithm as whuber suggested.

Using MATLAB

xs_name = "xs.txt";
ys_name = "ys.txt";

xs = arrayfun(@(x) str2double(x), xs);
xs = xs(~isnan(xs));

ys = arrayfun(@(y) str2double(y), ys);
ys = ys(~isnan(ys));

M = [xs'; ys']';

reduced_shape = reducepoly(M, 0.0001);

figure;
plot(xs, ys, 'r');
hold on;
plot(reduced_shape(:,1), reduced_shape(:, 2), 'b');


Original polygon in red.

Reduced polygon in blue.