1

I have a black-white image where I want to fit 2 separate lines to the edges in the image. It is easy to fit a single line with opencv, with the code below. How do I fit two best possible lines to this image. Here is the input image and 1 line result. I need something that is like the last image below.

def fit_line_to_edges(edges):
    points = np.where(edges>0)
    points = np.array([points[1], points[0]]).T
    [vx,vy,x,y] = cv2.fitLine(points, cv2.DIST_L2,0,0.01,0.01)
    lefty = int((-x*vy/vx) + y)
    righty = int(((edges.shape[1]-x)*vy/vx)+y)
    cv2.line(edges, (edges.shape[1]-1,righty), (0,lefty), color = (136,155,112), thickness = 3)
    return edges

enter image description here

enter image description here

Ozcan
  • 726
  • 4
  • 13
  • Your choice of these two lines seems completely arbitrary. What is the rationale ? If there is none, just split the array at a random point and perform two fits. –  Sep 26 '22 at 10:20
  • @YvesDaoust the logic is keep the total fit error of 2 lines at minimum, best 2 fits. – Ozcan Sep 26 '22 at 11:02
  • The keywords "piecewise linear regression" should bring a few useful results. – Stef Sep 26 '22 at 11:08
  • Also, there are variations of decision trees that end with linear regression instead of classification; and typically, these decision/regression tree algorithms can take the number of leaves as a parameter. See for instance this related question: [Decision tree with final decision being a linear regression](https://datascience.stackexchange.com/questions/65585/decision-tree-with-final-decision-being-a-linear-regression) – Stef Sep 26 '22 at 11:09
  • Also [Wikipedia: Segmented regression](https://en.wikipedia.org/wiki/Segmented_regression) – Stef Sep 26 '22 at 11:11
  • 1
    Also, if you call f(x) the regression error when using x as the cutting-point for the two pieces, then I expect f should be a unimodal function, that is, a function with a single global minimum, which is decreasing on the left of its global minimum and increasing on the right of its global minimum. You could use [golden-section search](https://en.wikipedia.org/wiki/Golden-section_search) to find the global minimum, ie the optimal cutting-point. (It's not completely obvious to me whether f should always be unimodal no matter your data, but at least for not too-unnice data, it should be mostly) – Stef Sep 26 '22 at 11:13
  • 2
    This related question defines a simple `piecewise_linear` custom function, then uses it with `scipy.optimise.curve_fit` to fit the data: [Curvefitting optimization error when fitting piecewise linear function](https://stackoverflow.com/questions/45091604/curvefitting-optimization-error-when-fitting-piecewise-linear-function). This is a bit brutal as `curve_fit` doesn't take advantage of the specific problem and uses a general optimisation algorithm, but it works. – Stef Sep 26 '22 at 11:21
  • Other similar questions: [How to apply piecewise linear fit in python?](https://stackoverflow.com/questions/29382903/how-to-apply-piecewise-linear-fit-in-python); [Piecewise linear fit with n breakpoints](https://stackoverflow.com/questions/46218934/piecewise-linear-fit-with-n-breakpoints); [Piecewise linear fit](https://stackoverflow.com/questions/61217320/piecewise-linear-fit); [Optimising a piecewise linear regression](https://stackoverflow.com/questions/56341347/optimizing-a-piecewise-linear-regression) – Stef Sep 26 '22 at 11:24
  • Anyway, read the doc for `scipy.optimise.curve_fit`: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html – Stef Sep 26 '22 at 11:25
  • This page should also prove interesting: [xavierdupre.fr : Piecewise linear regression with scikit-learn predictors](http://www.xavierdupre.fr/app/mlinsights/helpsphinx/notebooks/piecewise_linear_regression.html) – Stef Sep 26 '22 at 11:30
  • 1
    This question also discusses a variety of options and provides useful links: [segmented linear regression in python?](https://stackoverflow.com/questions/8999389/segmented-linear-regression-in-python) – Stef Sep 26 '22 at 11:33
  • I was thinking something similar to RANSAC, where select random points fit a line, then fit a second line to the rest of the points, repeat this many times and select the one with minimum total error – Ozcan Sep 26 '22 at 11:51
  • @Ozcan Yes that's a possibility. And it should work. Note that it doesn't make sense to select a completely random subset of points; rather, you should select a random cutting-point and take all the points to the left (or right). However, given my "guess" that the total error should be a unimodal function of the cutting-point, I think that searching for the cutting-point randomly is a bit wasteful. That's why I suggested [golden-section search](https://en.wikipedia.org/wiki/Golden-section_search) – Stef Sep 26 '22 at 12:01
  • @Stef Hi there. Do we need to use ordered points to use the piecewise linear function? I mean is it like use the first N points from left to right and fit a line, then use the rest of points and fit a second line? – Ozcan Sep 26 '22 at 12:02
  • Better perform a parabolic fit. :-) –  Sep 26 '22 at 12:43
  • Reminder: `log( exp( a ( x - c) )^p + exp( b ( x - c ) )^p ) / p + d` with `p > 0`is a nice continuous approximation for a piecewise linear function with slopes `a` and `b`, meeting at `( c, d )` – mikuszefski Sep 27 '22 at 05:46

1 Answers1

0

Thanks for the help in the comments. Here is the working code using a piecewise linear function

from scipy import optimize
import matplotlib.pyplot as plt
import numpy as np
import cv2

def fit_piecewise_lines(edges):
    def piecewise_linear(x, x0, y0, k1, k2):
        return np.piecewise(x, [x < x0], [lambda x:k1*x + y0-k1*x0, lambda x:k2*x + y0-k2*x0])

    # get x and y points from edges
    points = np.where(edges>0)
    points = np.array([points[1], points[0]]).T
    x = points[:,0]
    y = points[:,1]

    # initial point
    p_init = [np.mean(x), np.mean(y), 1, 1]

    # set initial parameter estimates
    p, e = optimize.curve_fit(piecewise_linear, x, y, p0=p_init)

    # get the point locations
    pt1 = (np.min(x).astype(int), piecewise_linear(np.min(x), *p).astype(int))
    pt3 = (np.max(x).astype(int), piecewise_linear(np.max(x), *p).astype(int))
    pt2 = (p[0].astype(int), p[1].astype(int))

    # plot it
    cv2.line(edges, pt1, pt2, (255,0,0), 5)
    cv2.line(edges, pt2, pt3, (0,255,0), 5)


image_file = "./edges.png"
edges = cv2.imread(image_file)

fit_piecewise_lines(edges)

plt.figure(figsize=(10,10))
plt.imshow(edges)

enter image description here

Ozcan
  • 726
  • 4
  • 13