How to estimate piece-wise smooth fit to a noisy mask?

Question

I have a binary mask that originated from a piece-wise smooth curve. Each underlying curve segment is smooth, but connections between segments may or may not be smooth. The mask is noisy so it might contain several pixels around the underlying curve. Image below shows an example of such input:

I want to estimate a fit to the underlying curve, given this input without any other prior knowledge, that will be able to provide both smooth and non-smooth connections.

I'm working in python so I tried several methods available there, such as numpy's polynomial fit, scipy's spline smoothing and pyqt-fit's non-parameteric regression, but was not able to get to the desired output. Here's a code example:

    from imageio import imread
    import random
    import numpy as np
    from scipy.interpolate import UnivariateSpline
    import pyqt_fit.nonparam_regression as smooth
    from pyqt_fit import npr_methods
    from matplotlib import pyplot as plt

    mask = imread(r'C:\mask.bmp')
    ys, xs = np.where(mask)
    # add tiny noise and sort - silly way to comply to UnivariateSpline's requirement of "x must be strictly increasing"
    xs = np.array([x + random.random() * 1e-4 for x in xs])
    sorter = np.argsort(xs)
    xs = xs[sorter]
    ys = ys[sorter]
    # polynomial fit
    p = np.poly1d(np.polyfit(xs, ys, 5))
    # spline smoothing
    spl = UnivariateSpline(xs, ys, k=3, s=1e9)
    # non-parameteric regression
    k = smooth.NonParamRegression(xs, ys, method=npr_methods.LocalPolynomialKernel(q=3))
    k.fit()

    plt.figure()
    plt.imshow(mask, cmap='gray')
    linexs = np.array(range(mask.shape[1]))
    plt.plot(linexs, k(linexs), 'y', lw=1)
    plt.plot(linexs, spl(linexs), 'g', lw=1)
    plt.plot(linexs, p(linexs), 'b', lw=1)
    plt.show()

For the parameters shown in this example, these fits both fail to capture the non-smooth connection on the left, as well as to provide a good fit for the "tail" on the right:

Expected results should behave like the red curve in the image below, where I expect the curve at location 1 to be non-smooth and at location 2 to be smooth.

I'd be happy to get a reference to a suitable algorithm. If there's also a python implementation, that would be a plus.

Since you are asking for an _algorithm_ it seems more like slightly theoretical computer science question more suitable for [Computer Science](https://cs.stackexchange.com/) community. — sophros, May 31 '19 at 12:43
Any solution of the form "y = f(x)" will yield a single value of y for any single value of x, so the "tail" on the right-hand side - which appears to require two different values of y for the same x - cannot be of this form. — James Phillips, May 31 '19 at 14:23
@JamesPhillips Right, I'm looking for a different kind of solution. — anemone, May 31 '19 at 19:44
Have a look at this question: https://stackoverflow.com/questions/52014197/how-to-interpolate-a-2d-curve-in-python — Christian K., Jun 05 '19 at 15:20
also probably related: https://stackoverflow.com/questions/3983613/find-tunnel-center-line — mikuszefski, Jun 13 '19 at 07:21

score 1 · Answer 1 · answered May 14 '20 at 06:52

This could be close to what you want. The trick is to first skeletonize the mask (see https://scikit-image.org/docs/dev/auto_examples/edges/plot_skeleton.html) to extract the points along which the curve to be estimated.

from skimage.morphology import skeletonize
import matplotlib.pyplot as plt
import cv2
import numpy as np
from scipy.interpolate import interp1d

img = cv2.imread("./data/bAiew.png", 0)
img = cv2.medianBlur(img, 11)
img = img//255

# Sparse skeleton
skeleton = skeletonize(img, method='lee')

# Into dense
contours, _ = cv2.findContours(skeleton, 0,cv2.CHAIN_APPROX_NONE)
arch = contours[0]
x,y = arch[...,0].squeeze(), arch[...,1].squeeze()

# Fitting a curve
xx, yy = x[0::15], y[0::15] #<- sample every 15th element to see that the interpolate really works
f = interp1d(xx, yy)

plt.figure(figsize=(16,9))
plt.subplot(221)
plt.imshow(img)
plt.title('original mask')

plt.subplot(222)
plt.imshow(skeleton)
plt.title('skeleton')

plt.subplot(223)
plt.scatter(xx,yy)
plt.ylim([img.shape[0],0])
plt.xlim([0, img.shape[1]])
plt.title('sampled points')

plt.subplot(224)
plt.plot(xx,f(xx))
plt.ylim([img.shape[0],0])
plt.xlim([0, img.shape[1]])
plt.title('interpolate')
plt.show()

Output:

How to estimate piece-wise smooth fit to a noisy mask?

1 Answers1