Блог

Допустим, у меня есть следующие данные (измерения):

Как вы можете видеть, есть много острых точек (т. Е. Там, где наклон сильно меняется). Поэтому было бы неплохо провести еще несколько измерений вокруг этих точек. Для этого я написал скрипт:

1. Я вычисляю кривизну 3 последовательных точек: Кривизна Менгера:https://en.wikipedia.org/wiki/Menger_curvature#Definition

2. Затем я решаю, какие значения я должен выполнить повторную выборку, основываясь на кривизне.

…и я повторяю, пока средняя кривизна не уменьшится… но это не работает, потому что оно увеличивается. Вы знаете почему?

Вот полный код (остановлен после того, как длина значений x достигнет 60):

 import numpy as np
import matplotlib.pyplot as plt

def curvature(A,B,C):
    """Calculates the Menger curvature fro three Points, given as numpy arrays.
    Sources:
    Menger curvature: https://en.wikipedia.org/wiki/Menger_curvature#Definition
    Area of a triangle given 3 points: https://math.stackexchange.com/questions/516219/finding-out-the-area-of-a-triangle-if-the-coordinates-of-the-three-vertices-are
    """

    # Pre-check: Making sure that the input points are all numpy arrays
    if any(x is not np.ndarray for x in [type(A),type(B),type(C)]):
        print("The input points need to be a numpy array, currently it is a ", type(A))

    # Augment Columns
    A_aug = np.append(A,1)
    B_aug = np.append(B,1)
    C_aug = np.append(C,1)

    # Caclulate Area of Triangle
    matrix = np.column_stack((A_aug,B_aug,C_aug))
    area = 1/2*np.linalg.det(matrix)

    # Special case: Two or more points are equal 
    if np.all(A == B) or  np.all(B == C):
        curvature = 0
    else:
        curvature = 4*area/(np.linalg.norm(A-B)*np.linalg.norm(B-C)*np.linalg.norm(C-A))

    # Return Menger curvature
    return curvature


def values_to_calulate(x,curvature_list, max_curvature):
    """Calculates the new x values which need to be calculated
    Middle point between the three points that were used to calculate the curvature """
    i = 0
    new_x = np.empty(0)
    for curvature in curvature_list:
        if curvature > max_curvature:
            new_x = np.append(new_x, x[i] (x[i 2]-x[i])/3 )
        i = i 1
    return new_x


def plot(x,y, title, xLabel, yLabel):
    """Just to visualize"""

    # Plot
    plt.scatter(x,y)
    plt.plot(x, y, '-o')

    # Give a title for the sine wave plot
    plt.title(title)

    # Give x axis label for the sine wave plot
    plt.xlabel(xLabel)

    # Give y axis label for the sine wave plot
    plt.ylabel(yLabel)
    plt.grid(True, which='both')
    plt.axhline(y=0, color='k')


    # Display the sine wave
    plt.show
    plt.pause(0.05)

### STARTS HERE


# Get x values of the sine wave
x = np.arange(0, 10, 1);

# Amplitude of the sine wave is sine of a variable like time
def function(x):
    return 1 np.sin(x)*np.cos(x)**2
y = function(x)

# Plot it
plot(x,y, title='Data', xLabel='Time', yLabel='Amplitude')


continue_Loop = True

while continue_Loop == True :
    curvature_list = np.empty(0)
    for i in range(len(x)-2):
        # Get the three points
        A = np.array([x[i],y[i]])
        B = np.array([x[i 1],y[i 1]])
        C = np.array([x[i 2],y[i 2]])

        # Calculate the curvature
        curvature_value = abs(curvature(A,B,C))
        curvature_list = np.append(curvature_list, curvature_value)



    print("len: ", len(x) )
    print("average curvature: ", np.average(curvature_list))

    # Calculate the points that need to be added 
    x_new = values_to_calulate(x,curvature_list, max_curvature=0.3)

    # Add those values to the current x list:
    x = np.sort(np.append(x, x_new))

    # STOPED IT AFTER len(x) == 60
    if len(x) >= 60:
        continue_Loop = False

    # Amplitude of the sine wave is sine of a variable like time
    y = function(x)

    # Plot it
    plot(x,y, title='Data', xLabel='Time', yLabel='Amplitude')
  

Вот как это должно выглядеть:

Редактировать:

Если вы позволите ему работать еще дальше… :

Итак, подведите итог моим комментариям выше:

• вы вычисляете среднюю кривизну вашей кривой, у которой нет причин переходить к 0. В каждой точке, независимо от того, насколько близки ваши точки, радиус окружности будет сходиться к любой кривизне в этой точке, а не к 0.

• альтернативой было бы использовать изменение абсолютной производной между двумя точками: продолжайте выборку до abs(d(df/dx)) < some_threshold того момента, когда d(df/dx) = (df/dx)[n] - (df/dx)[n-1]