Manera numérica de resolver la curvatura de una curva.

Question

Soy consciente de la ecuación para resolver la curvatura de una curva. Es$\kappa = \displaystyle\frac{|r' \times r''|}{|r'|^3}$, donde$r$ es la parametrización de la curva.

En ausencia de la ecuación paramétrica, donde solo se dan los puntos de datos, ¿cómo puede resolver / aproximar la curvatura con un error sobre$O(h^2)$?

Answer 1

Como dbx se menciona en los comentarios, puede utilizar el radio de un círculo que pasa por tres puntos. La curvatura es la inversa de la radio. La radio puede ser encontrado a través de $$ r = \frac{abc}{4k}, $$ donde $a$, $b$ y $c$ indican las distancias entre los tres puntos y $k$ denota el área del triángulo formado por los tres puntos. Obviamente, la curvatura es el recíproco de este, por lo tanto $$ \kappa = \frac{4k}{abc} $$ Se me ocurrió este código en el pasado en python. A continuación está el código. He utilizado la fórmula de Herón para calcular el área del triángulo $k$.

# points_utm is a 3-by-2 array, containing the easting and northing coordinates of 3 points
# Compute distance to each point
a = np.hypot(points_utm[0, 0] - points_utm[1, 0], points_utm[0, 1] - points_utm[1, 1])
b = np.hypot(points_utm[1, 0] - points_utm[2, 0], points_utm[1, 1] - points_utm[2, 1])
c = np.hypot(points_utm[2, 0] - points_utm[0, 0], points_utm[2, 1] - points_utm[0, 1])

# Compute inverse radius of circle using surface of triangle (for which Heron's formula is used)
k = np.sqrt((a+(b+c))*(c-(a-b))*(c+(a-b))*(a+(b-c))) / 4    # Heron's formula for triangle's surface
den = a*b*c  # Denumerator; make sure there is no division by zero.
if den == 0.0:  # Very unlikely, but just to be sure
    curvature = 0.0
else:
    curvature = 4*k / den

Nota adicional: este podría no ser óptima ya que sólo tres puntos se utilizan para calcular la curvatura. Si más puntos están disponibles, sería beneficioso para el uso de ellos. Sin embargo, se puede utilizar el mismo enfoque: ajuste de un círculo a través de los puntos de datos y calcular el recíproco de la radio. Sin embargo, no es trivial para encontrar el mejor círculo de ajuste para los datos. Para obtener más información, consulte http://scipy-cookbook.readthedocs.io/items/Least_Squares_Circle.html.

Basado en el anterior enlace, escribí un pequeño script en python que se muestra cómo se puede calcular la curvatura utilizando un lineal de mínimos cuadrados. Aquí está el código:

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


class ComputeCurvature:
    def __init__(self):
        """ Initialize some variables """
        self.xc = 0  # X-coordinate of circle center
        self.yc = 0  # Y-coordinate of circle center
        self.r = 0   # Radius of the circle
        self.xx = np.array([])  # Data points
        self.yy = np.array([])  # Data points

    def calc_r(self, xc, yc):
        """ calculate the distance of each 2D points from the center (xc, yc) """
        return np.sqrt((self.xx-xc)**2 + (self.yy-yc)**2)

    def f(self, c):
        """ calculate the algebraic distance between the data points and the mean circle centered at c=(xc, yc) """
        ri = self.calc_r(*c)
        return ri - ri.mean()

    def df(self, c):
        """ Jacobian of f_2b
        The axis corresponding to derivatives must be coherent with the col_deriv option of leastsq"""
        xc, yc = c
        df_dc = np.empty((len(c), x.size))

        ri = self.calc_r(xc, yc)
        df_dc[0] = (xc - x)/ri                   # dR/dxc
        df_dc[1] = (yc - y)/ri                   # dR/dyc
        df_dc = df_dc - df_dc.mean(axis=1)[:, np.newaxis]
        return df_dc

    def fit(self, xx, yy):
        self.xx = xx
        self.yy = yy
        center_estimate = np.r_[np.mean(xx), np.mean(yy)]
        center = optimize.leastsq(self.f, center_estimate, Dfun=self.df, col_deriv=True)[0]

        self.xc, self.yc = center
        ri = self.calc_r(*center)
        self.r = ri.mean()

        return 1 / self.r  # Return the curvature


# Apply code for an example
x = np.r_[36, 36, 19, 18, 33, 26]
y = np.r_[14, 10, 28, 31, 18, 26]
comp_curv = ComputeCurvature()
curvature = comp_curv.fit(x, y)

# Plot the result
theta_fit = np.linspace(-np.pi, np.pi, 180)
x_fit = comp_curv.xc + comp_curv.r*np.cos(theta_fit)
y_fit = comp_curv.yc + comp_curv.r*np.sin(theta_fit)
plt.plot(x_fit, y_fit, 'k--', label='fit', lw=2)
plt.plot(x, y, 'ro', label='data', ms=8, mec='b', mew=1)
plt.xlabel('x')
plt.ylabel('y')
plt.title('curvature = {:.3e}'.format(curvature))
plt.show()