# Approximating the boundary between two sets of points (in 2D): Fitting a region

Given two sets of points $$p_{\text{in},i}$$ and $$p_{\text{out},j}$$ inside and outside of what I intuitively call a "region", I would like to estimate and describe the boundary of this region. Two kinds of answer are welcome:

• A parametrized description of the boundary
• A estimation function that predicts whether a new point is inside the region.

So far, I have found lots of clustering algorithms, which may be related, but I do not see how to apply them here. Here, in/out-ness is part of the input, not the result.

I arrived at this question from an experimental context. A nice follow-up question is "where to probe next" (i.e. at which new location would in/out-knowledge change the model most? In the image below: bottom left)

My first naive approach is something like "a circle that minimizes some kind of signed distance to all points", but I want more complex shapes.

Some more thoughts:

• I'm solely interested in 2D right now
• I'm fairly proficient at optimization, even pointers to suitable parametrizable "descriptions of regions" would be welcome. I know these two sequences, for example, but they don't seem suitably bounded.
• The number of optimization parameters should be scalable: for just one point in and one point inside, even returning a circle is overfitting, but for the data below, a more complex shape of the region makes sense.
• There will be outliers
• An alogrithm for convex regions only would be acceptable, but a more general one is preferred.
• Weights (uncertainties) for the points are available.

My original question was unclear in some ways. Thanks to the first answers, I can now specify some things more clearly:

• "What I intuitively call a region" has two important properties: it is finite and connected
• There will be outliers: Some points will be wrongly classified, and the model shall not overfit to include all points. Hence "number of parameters should be scalable" to avoid/adjust overfitting, just like one will usually not fit a polynomial of degree $$n$$ to $$n-1$$ points.

• If you can provide some data, I can test my idea and see if it matches what you want Nov 21 at 5:39

# Theory

Clustering is unlikely to work in this case because your red points are separated from each other by the green points. You could use more clusters, but this will require a lot of manual inspection and fiddling.

A standard approach to this sort of problem is to use a nonlinear support vector machine. The idea, simply described, is that although your data may not be separable in two dimensions projecting the data into a higher dimension always guarantees a clean separator which can then be used to classify future points of interest. Furthermore, we can choose this separator so that the distance between the separator and any of the data points is maximized.

Here are a couple of images showing the idea:

# Results

I used WebPlotDigitizer to extract your data in arbitrary, but aspect-ratio-preserving units and then ran this through a support-vector machine using a radial basis function kernel. This gives the following separator:

Notice that, as promised, this cleanly separates your data into two classes. The separator itself is the solid grey line while the margins around the separator (the distance between the separator and the data) are shown with dotted lines. The distance between the separator and these margins is the maximum versus all other choices for this kernel. The data points which induce the margins are circled.

If we add a hundred new points and classify them using our SVM, we get this:

# Implementation

I use Python to implement the above ideas.

import matplotlib.pyplot as plt
import numpy as np
from sklearn import svm

data = [
[ 0.7490599550363091 ,  9.24443743612264   , 0],
[ 4.12034765121126   ,  8.885539916739752  , 0],
[ 7.6539862433328665 ,  9.612037680925432  , 0],
[ 7.914094745271155  ,  0.9253511681985884 , 0],
[ 8.994538029638242  ,  0.5077641947442704 , 0],
[ 9.371974940345583  ,  9.608153819994236  , 0],
[12.320128383367402  ,  6.088728506174032  , 0],
[12.411426660363531  ,  4.715352126892222  , 0],
[12.558352293037347  , 11.229593835418646  , 0],
[13.114336626127825  ,  5.703722810531428  , 0],
[13.381386496538893  ,  9.854563219073547  , 0],
[13.886274645038084  ,  8.480251836234228  , 0],
[14.13510342320811   ,  6.595572357695316  , 0],
[14.463275899337807  ,  9.915985760466931  , 0],
[15.05984795641495   ,  7.296033868971857  , 0],
[ 3.096302985685499  ,  6.907935469254409  , 1],
[ 4.878692156862327  ,  7.159379502874164  , 1],
[ 7.4911533079089345 ,  8.366972559074298  , 1],
[ 7.986689890548966  ,  3.89506631355318   , 1],
[ 8.498712223311847  ,  4.883868537295852  , 1],
[ 9.869150457208352  ,  5.679125024633843  , 1],
[10.527327159481406  ,  2.3884159656508306 , 1],
[10.671071331605198  ,  7.848836741512311  , 1],
[12.788477939489049  ,  3.1497246315161327 , 1],
[12.86531503216689   ,  7.524534353757313  , 1],
]

# Load data into numpy
data = np.array(data)

# Separate data into x and y values; predictors and observations
data_x = data[:,0:2]  # x-value/predictor
data_y = data[:,2]    # y-value/observation

# In these "learning" style problems you often want to avoid over-fitting, so
# it might make sense to split your data into training and test sets. Here we
# simply use all data as the training data.
x_train = data_x[:,:]
y_train = data_y[:]

# Using the RBF (radial basis function) kernel with an SVM projects the points
# into a higher-dimensional space where a clean boundary is possible and then
# lowers that boundary back into the space the points live in.
model = svm.SVC(kernel='rbf', C=1e6)
model.fit(x_train, y_train)

def plot_model(data_x, data_y, predict_x=None, predict_y=None) -> None:
# Let's plot some stuff
fig, ax = plt.subplots(figsize=(12, 7))

# Create grid to evaluate model
xx = np.linspace(-1, max(data_x[:,0]) + 1, len(data_x))
yy = np.linspace(0, max(data_x[:,1]) + 1, len(data_x))
YY, XX = np.meshgrid(yy, xx)
xy = np.vstack([XX.ravel(), YY.ravel()]).T

# Assigning different colors to the classes
colors = data_y
colors = np.where(colors == 1, '#8C7298', '#4786D1')
ax.scatter(data_x[:,0], data_x[:,1],c=colors)

if predict_x is not None:
assert predict_y is not None
colors = predict_y
colors = np.where(colors == 1, '#8C7298', '#4786D1')
ax.scatter(predict_x[:,0], predict_x[:,1],c=colors,s=4)

# Get the separator
Z = model.decision_function(xy).reshape(XX.shape)

# Draw the decision boundary and margins
ax.contour(XX, YY, Z, colors='k', levels=[-1, 0, 1], alpha=0.5, linestyles=['--', '-', '--'])

# Highlight support vectors with a circle around them
ax.scatter(model.support_vectors_[:, 0], model.support_vectors_[:, 1], s=100, linewidth=1, facecolors='none', edgecolors='k')

plot_model(data_x, data_y)
plt.show()

# Generate some previously-unseen data
new_x = np.random.uniform(low=min(data_x[:,0]), high=max(data_x[:,0]), size=(100,1))
new_y = np.random.uniform(low=min(data_x[:,1]), high=max(data_x[:,1]), size=(100,1))
new_data = np.hstack((new_x, new_y))
predictions_for_new_data = model.predict(new_data)

plot_model(data_x, data_y, new_data, predictions_for_new_data)
plt.show()

• Amazing! Much better than my answer, and I learned something today. Thanks :) Nov 24 at 5:26
• Thanks for this well-written answer! Although I would be more comfortable with an analytic set of base functions, this looks like a promising approach. However, the results are very overfitted: A large number of points are exactly on the 1/-1 isolines. Given that the data includes outliers, some "training" points should end up on the wrong side of the classifier. I'm only superficially familiar with SVMs - is there a way to reduce the number of DOF to get a "close enough, but less complex" answer? Also, is there a way to force connectedness of the result? Nov 25 at 7:15
• @Lichtkarussell: You should expect a large number of points to be on the 1/-1 isolines, that's why the lines are drawn there - they are maximizing the margin and the maximal margin must intersect points. "Close enough, but less complex" is a long way of saying regularization. The docs indicate a parameter C that can be adjusted for this. I don't know if you can force connectedness. In general, probably not, since continuous high-dimensional manifolds may not be continuous in lower dimensions. Nov 25 at 7:45
• Thanks! Found C in the docs, gives a pretty good result. There's a fine point between "on the isolines" and "near the isolines" - your plots looked like the former, but are in fact the latter, with C=1e6 yielding subplixel-nearness. With C=1 (and after realzing the SVM only works if I shift my points to the origin), I get pretty decent-looking results. Please forgive me for not choosing this as an answer (yet?) - I'd really like something that is finite and connected by definition, not just by chance for some inputs. Nov 25 at 8:14

I don't know if the following is a good idea. But it is an idea, and I hope that it helps.

This problem can be recast to "find a function $$f: \mathbb{R}^2\to \mathbb{R}$$ and $$z\in\mathbb{R}$$ s.t. $$f(x,y)-z \geq 0$$ if the point $$(x,y)$$ is "inside", and negative otherwise. For the basis of $$f$$ you can pick tensor products of Legendre polynomials up to degree 10 as an example. At this point, you have an LP problem to determine the coefficients of the basis functions and the constant $$z$$, which you can solve with many solvers. Then the prediction is as simple evaluating $$f$$ at a point and comparing to $$z$$.

Using an idea from Richard (see their answer at https://scicomp.stackexchange.com/a/40497/33791 ), I extracted the points from the graph you provided. I assumed the bottom left point is $$(0,0)$$ and the top right point was $$(1,1)$$. If not, consider that I worked with normalized data.

I implemented my idea in MATLAB, if of interest I can clean up and share the code. This is the result:

The surface I obtained may have some undesirable properties. For example, the data does not indicate that top right and top left parts would be interior. But my solution predicts that those areas would be in the interior of the region.

On the other hand, this method can clearly generate very complex, non-convex and disconnected regions, which is nice. It can also deal with uncertainties. For example, in my model:

• if $$f(x,y)\ge\epsilon>0$$ then $$(x,y)$$ is inside,
• if $$\epsilon>f(x,y)>0$$ then $$(x,y)$$ is likely to be inside,
• if $$0>f(x,y)>-\epsilon$$ then $$(x,y)$$ is likely to be outside, and
• if $$0>-\epsilon\ge f(x,y)$$ then $$(x,y)$$ is outside.

You can see this in the artifacts of the plot, as the boundary is very sharp. But by setting the tolerance larger than I did, it should be possible to have a softer gradient, hence more blurred boundary.

• Yes, the "find a function" approach seems like a good idea (e.g. f(x,y)=x²+y²-c for a cricular boundary). However, Legendre polynomials are a poor choice because they will always have oscillating values at the borders (I think?), so they cannot describe a connected region. This shows in your implementation, too. (Along with the problem of arbitrarily choosing the bounds - it may well be that the "true" region extends below the bottom left in the data) The Zernike polynomials are linked in the Q - they are more finite, but oscillate, too. With a suitable base, this might work, though. Nov 25 at 7:12