.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/plot_splinebox_vs_scipy_line.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_plot_splinebox_vs_scipy_line.py: Comparison splinebox and scipy: edge fitting -------------------------------------------- This example compares splinebox and scipy when trying to fit a spline to an end of an image. .. GENERATED FROM PYTHON SOURCE LINES 8-20 .. code-block:: Python # sphinx_gallery_thumbnail_number = 2 import cycler import matplotlib as mpl import matplotlib.pyplot as plt import numpy as np import scipy import scipy.optimize import skimage import splinebox .. GENERATED FROM PYTHON SOURCE LINES 21-24 First, we will set some defaults to style our plot. You can ingnore this section if you don't care about the style of the plots. .. GENERATED FROM PYTHON SOURCE LINES 24-28 .. code-block:: Python mpl.rcParams["image.cmap"] = "gray" mpl.rcParams["lines.linewidth"] = 2 mpl.rcParams["axes.prop_cycle"] = cycler.cycler(color=("r",)) .. GENERATED FROM PYTHON SOURCE LINES 29-30 In this example, we will use a crop of the scikit-image's `text` example image. .. GENERATED FROM PYTHON SOURCE LINES 30-35 .. code-block:: Python img = skimage.data.text() img = img[45:78, 300:430] plt.imshow(img) plt.show() .. image-sg:: /auto_examples/images/sphx_glr_plot_splinebox_vs_scipy_line_001.png :alt: plot splinebox vs scipy line :srcset: /auto_examples/images/sphx_glr_plot_splinebox_vs_scipy_line_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 36-40 Our goal is to fit a spline to the integral sign on the page. We won't cover how to automatically find a good initial guess but instead just select five pixels that are more or less equally spaced along the integral as the initial not of our spline. .. GENERATED FROM PYTHON SOURCE LINES 40-42 .. code-block:: Python initial_knots = np.array([[24, 12], [23, 47], [16, 72], [8, 97], [8, 124]]) .. GENERATED FROM PYTHON SOURCE LINES 43-48 In order to evaluate how well our spline fits the line, we will have to interpolate the pixel values at non-integer position. Since the line is black on a white background, our goal is to move the spline in a way that minimises the average pixel value under it. This is commonly refered to as the image energy. .. GENERATED FROM PYTHON SOURCE LINES 48-50 .. code-block:: Python interpolator = scipy.interpolate.RectBivariateSpline(np.arange(img.shape[0]), np.arange(img.shape[1]), img, s=1) .. GENERATED FROM PYTHON SOURCE LINES 51-53 We will start by constructing an initial spline using splinebox. .. GENERATED FROM PYTHON SOURCE LINES 53-58 .. code-block:: Python M = len(initial_knots) basis_function = splinebox.B3() spline = splinebox.Spline(M, basis_function, closed=False) spline.knots = initial_knots .. GENERATED FROM PYTHON SOURCE LINES 59-61 Let's define the parameter values at which we want to sample the spline. Here, we chose to sample 50 points inbetween two knots. .. GENERATED FROM PYTHON SOURCE LINES 61-63 .. code-block:: Python t = np.linspace(0, M - 1, M * 50) .. GENERATED FROM PYTHON SOURCE LINES 64-66 In order to compar the fitted spline to the intial one, we save it's positions and knots for plotting later on. .. GENERATED FROM PYTHON SOURCE LINES 66-70 .. code-block:: Python initial_vals = spline.eval(t) initial_knots = spline.knots .. GENERATED FROM PYTHON SOURCE LINES 71-80 Before we can start fitting the spline, we have to define the loss function. In this example our loss function consists of two parts, the image energy, discussed above and an internal energy. Here, we choose to use the curvilinear reparametrization energy as our internal energy. It promotes equidistant spacing of the knots in terms of arc length. In practice, this avoids sharp bends and stops the spline from looping/folding back on itself. Without it, the image energy would reward the spline for visiting the darkest pixels multiple times. The parameter alpha can be used balance the contribution of the image and internal energies. .. GENERATED FROM PYTHON SOURCE LINES 80-88 .. code-block:: Python def loss_function_splinebox(control_points, alpha): spline.control_points = control_points.reshape((-1, 2)) coordinates = spline.eval(t) image_energy = np.mean(interpolator(coordinates[:, 0], coordinates[:, 1], grid=False)) internal_energy = spline.curvilinear_reparametrization_energy() return image_energy + alpha * internal_energy .. GENERATED FROM PYTHON SOURCE LINES 89-91 For the fitting procedure we can simply use scipy. The value for alpha was empirically set to 500. .. GENERATED FROM PYTHON SOURCE LINES 91-94 .. code-block:: Python initial_control_points = spline.control_points scipy.optimize.minimize(loss_function_splinebox, initial_control_points.flatten(), args=(500,)) .. rst-class:: sphx-glr-script-out .. code-block:: none message: Optimization terminated successfully. success: True status: 0 fun: 36.84064228311857 x: [-1.546e+01 1.866e+01 ... 2.018e+01 5.141e+01] nit: 119 jac: [ 9.537e-07 4.768e-07 ... 0.000e+00 4.768e-07] hess_inv: [[ 1.282e+02 -6.134e+01 ... 1.102e+01 1.306e+01] [-6.134e+01 3.748e+02 ... -2.810e+01 8.961e+01] ... [ 1.102e+01 -2.810e+01 ... 1.693e+02 -9.179e+01] [ 1.306e+01 8.961e+01 ... -9.179e+01 3.397e+02]] nfev: 2055 njev: 137 .. GENERATED FROM PYTHON SOURCE LINES 95-96 Let's plot the results. .. GENERATED FROM PYTHON SOURCE LINES 96-112 .. code-block:: Python fitted_vals = spline.eval(t) fitted_knots = spline.knots fix, axes = plt.subplots(2, 1, sharex=True, sharey=True) axes[0].imshow(img) axes[0].plot(initial_vals[:, 1], initial_vals[:, 0], label="initial spline") axes[0].scatter(initial_knots[:, 1], initial_knots[:, 0], label="initial knots") axes[0].legend() axes[1].imshow(img) axes[1].plot(fitted_vals[:, 1], fitted_vals[:, 0], label="fitted spline") axes[1].scatter(fitted_knots[:, 1], fitted_knots[:, 0], label="fitted knots") axes[1].legend() plt.suptitle("SplineBox") plt.tight_layout() plt.show() .. image-sg:: /auto_examples/images/sphx_glr_plot_splinebox_vs_scipy_line_002.png :alt: SplineBox :srcset: /auto_examples/images/sphx_glr_plot_splinebox_vs_scipy_line_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 113-117 Next, we will try to acchive the same fit using scipy. As before, we begin by construcing an initial spline using the initial knots. NOTE: you have to set `bc_type` in order to get a spline with the desired number of knots. .. GENERATED FROM PYTHON SOURCE LINES 117-125 .. code-block:: Python k = 3 t_knots = np.arange(M) spline = scipy.interpolate.make_interp_spline(t_knots, initial_knots, k=3, bc_type="natural") initial_vals = spline(t) initial_knots = spline(spline.t)[k:-k] .. GENERATED FROM PYTHON SOURCE LINES 126-128 Since scipy does not provide a function to compute the curvilinear reparametrization energy, we have to do it ourselfs. .. GENERATED FROM PYTHON SOURCE LINES 128-144 .. code-block:: Python def loss_function_scipy(control_points, alpha): spline.c = control_points.reshape((-1, 2)) coordinates = spline(t) image_energy = np.mean(interpolator(coordinates[:, 0], coordinates[:, 1], grid=False)) derivative = spline.derivative() integral = scipy.integrate.quad(lambda t: np.linalg.norm(derivative(t)), 0, M - 1) length = integral[0] c = (length / M) ** 2 integral = scipy.integrate.quad(lambda t: (np.linalg.norm(derivative(t)) ** 2 - c) ** 2, 0, M - 1) internal_energy = integral[0] / length**4 return image_energy + alpha * internal_energy initial_control_points = spline.c scipy.optimize.minimize(loss_function_scipy, initial_control_points.flatten(), args=(500,)) .. rst-class:: sphx-glr-script-out .. code-block:: none message: Optimization terminated successfully. success: True status: 0 fun: 36.84064228312302 x: [ 2.409e+01 1.117e+01 ... 6.737e+00 1.239e+02] nit: 68 jac: [ 1.431e-06 0.000e+00 ... 1.907e-06 2.861e-06] hess_inv: [[ 9.638e-01 1.367e-01 ... 5.749e-02 -3.243e-03] [ 1.367e-01 9.017e-01 ... -5.137e-03 -1.333e-02] ... [ 5.749e-02 -5.137e-03 ... 1.569e+00 -4.038e-01] [-3.243e-03 -1.333e-02 ... -4.038e-01 1.240e+00]] nfev: 1215 njev: 81 .. GENERATED FROM PYTHON SOURCE LINES 145-146 Let's take a look at the results. .. GENERATED FROM PYTHON SOURCE LINES 146-161 .. code-block:: Python fitted_vals = spline(t) fitted_knots = spline(spline.t[k:-k]) fix, axes = plt.subplots(2, 1, sharex=True, sharey=True) axes[0].imshow(img, cmap="gray") axes[0].plot(initial_vals[:, 1], initial_vals[:, 0], label="initial spline") axes[0].scatter(initial_knots[:, 1], initial_knots[:, 0], label="initial knots") axes[0].legend() axes[1].imshow(img, cmap="gray") axes[1].plot(fitted_vals[:, 1], fitted_vals[:, 0], label="fitted spline") axes[1].scatter(fitted_knots[:, 1], fitted_knots[:, 0], label="fitted knots") axes[1].legend() plt.suptitle("SciPy") plt.tight_layout() plt.show() .. image-sg:: /auto_examples/images/sphx_glr_plot_splinebox_vs_scipy_line_003.png :alt: SciPy :srcset: /auto_examples/images/sphx_glr_plot_splinebox_vs_scipy_line_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 34.547 seconds) .. _sphx_glr_download_auto_examples_plot_splinebox_vs_scipy_line.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_splinebox_vs_scipy_line.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_splinebox_vs_scipy_line.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_splinebox_vs_scipy_line.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_