Finding the corners of noisy polygons

Question

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.

I read about Lasso and decided to try that.

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

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

xs = np.loadtxt(xs_name).reshape(-1, 1)
ys = np.loadtxt(ys_name).reshape(-1, 1)

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")
ys <- read.table('ys.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?

Answer 1

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.)

Answer 2

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 = readlines(xs_name);
xs = arrayfun(@(x) str2double(x), xs);
xs = xs(~isnan(xs));

ys = readlines(ys_name);
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.