Алгоритм нахождения локальных максимумов

Question

У меня есть данные, которые всегда выглядят примерно так:

альтернативный текст http://michaelfogleman.com/static/images/chart.png

Мне нужен алгоритм для определения трех пиков.

Ось X на самом деле является положением камеры, а ось Y - это мера фокусировки / контраста изображения в этой позиции. Есть объекты на трех разных расстояниях, которые могут быть в фокусе, и мне нужно определить значения x для этих трех точек.

Средний горб всегда немного сложнее, даже для человека.

У меня есть самодельный алгоритм, который в основном работает, но мне интересно, есть ли стандартный способ извлечения локальных максимумов из функции, в которой может быть немного шума? Хотя пики легко преодолевают шум.

Кроме того, для данных камеры может быть полезен алгоритм, не требующий сканирования полного диапазона.

Редактировать : публикация кода Python, который я в итоге использовал. Он использует мой оригинальный код, который находит максимумы при заданном пороге поиска, и выполняет двоичный поиск, чтобы найти порог, который приводит к требуемому количеству максимумов.

Редактировать : Пример данных, включенных в код ниже. Новый код O (n) вместо O (n ^ 2).

def find_n_maxima(data, count):
    low = 0
    high = max(data) - min(data)
    for iteration in xrange(100): # max iterations
        mid = low + (high - low) / 2.0
        maxima = find_maxima(data, mid)
        if len(maxima) == count:
            return maxima
        elif len(maxima) < count: # threshold too high
            high = mid
        else: # threshold too low
            low = mid
    return None # failed

def find_maxima(data, threshold):
    def search(data, threshold, index, forward):
        max_index = index
        max_value = data[index]
        if forward:
            path = xrange(index + 1, len(data))
        else:
            path = xrange(index - 1, -1, -1)
        for i in path:
            if data[i] > max_value:
                max_index = i
                max_value = data[i]
            elif max_value - data[i] > threshold:
                break
        return max_index, i
    # forward pass
    forward = set()
    index = 0
    while index < len(data) - 1:
        maximum, index = search(data, threshold, index, True)
        forward.add(maximum)
        index += 1
    # reverse pass
    reverse = set()
    index = len(data) - 1
    while index > 0:
        maximum, index = search(data, threshold, index, False)
        reverse.add(maximum)
        index -= 1
    return sorted(forward & reverse)

data = [
    1263.900, 1271.968, 1276.151, 1282.254, 1287.156, 1296.513,
    1298.799, 1304.725, 1309.996, 1314.484, 1321.759, 1323.988,
    1331.923, 1336.100, 1340.007, 1340.548, 1343.124, 1353.717,
    1359.175, 1364.638, 1364.548, 1357.525, 1362.012, 1367.190,
    1367.852, 1376.275, 1374.726, 1374.260, 1392.284, 1382.035,
    1399.418, 1401.785, 1400.353, 1418.418, 1420.401, 1423.711,
    1425.214, 1436.231, 1431.356, 1435.665, 1445.239, 1438.701,
    1441.988, 1448.930, 1455.066, 1455.047, 1456.652, 1456.771,
    1459.191, 1473.207, 1465.788, 1488.785, 1491.422, 1492.827,
    1498.112, 1498.855, 1505.426, 1514.587, 1512.174, 1525.244,
    1532.235, 1543.360, 1543.985, 1548.323, 1552.478, 1576.477,
    1589.333, 1610.769, 1623.852, 1634.618, 1662.585, 1704.127,
    1758.718, 1807.490, 1852.097, 1969.540, 2243.820, 2354.224,
    2881.420, 2818.216, 2552.177, 2355.270, 2033.465, 1965.328,
    1824.853, 1831.997, 1779.384, 1764.789, 1704.507, 1683.615,
    1652.712, 1646.422, 1620.593, 1620.235, 1613.024, 1607.675,
    1604.015, 1574.567, 1587.718, 1584.822, 1588.432, 1593.377,
    1590.533, 1601.445, 1667.327, 1739.034, 1915.442, 2128.835,
    2147.193, 1970.836, 1755.509, 1653.258, 1613.284, 1558.576,
    1552.720, 1541.606, 1516.091, 1503.747, 1488.797, 1492.021,
    1466.720, 1457.120, 1462.485, 1451.347, 1453.224, 1440.477,
    1438.634, 1444.571, 1428.962, 1431.486, 1421.721, 1421.367,
    1403.461, 1415.482, 1405.318, 1399.041, 1399.306, 1390.486,
    1396.746, 1386.178, 1376.941, 1369.880, 1359.294, 1358.123,
    1353.398, 1345.121, 1338.808, 1330.982, 1324.264, 1322.147,
    1321.098, 1313.729, 1310.168, 1304.218, 1293.445, 1285.296,
    1281.882, 1280.444, 1274.795, 1271.765, 1266.857, 1260.161,
    1254.380, 1247.886, 1250.585, 1246.901, 1245.061, 1238.658,
    1235.497, 1231.393, 1226.241, 1223.136, 1218.232, 1219.658,
    1222.149, 1216.385, 1214.313, 1211.167, 1208.203, 1206.178,
    1206.139, 1202.020, 1205.854, 1206.720, 1204.005, 1205.308,
    1199.405, 1198.023, 1196.419, 1194.532, 1194.543, 1193.482,
    1197.279, 1196.998, 1194.489, 1189.537, 1188.338, 1184.860,
    1184.633, 1184.930, 1182.631, 1187.617, 1179.873, 1171.960,
    1170.831, 1167.442, 1177.138, 1166.485, 1164.465, 1161.374,
    1167.185, 1174.334, 1186.339, 1202.136, 1234.999, 1283.328,
    1347.111, 1679.050, 1927.083, 1860.902, 1602.791, 1350.454,
    1274.236, 1207.727, 1169.078, 1138.025, 1117.319, 1109.169,
    1080.018, 1073.837, 1059.876, 1050.209, 1050.859, 1035.003,
    1029.214, 1024.602, 1017.932, 1006.911, 1010.722, 1005.582,
    1000.332, 998.0721, 992.7311, 992.6507, 981.0430, 969.9936,
    972.8696, 967.9463, 970.1519, 957.1309, 959.6917, 958.0536,
    954.6357, 954.9951, 947.8299, 953.3991, 949.2725, 948.9012,
    939.8549, 940.1641, 942.9881, 938.4526, 937.9550, 929.6279,
    935.5402, 921.5773, 933.6365, 918.7065, 922.5849, 939.6088,
    911.3251, 923.7205, 924.8227, 911.3192, 936.7066, 915.2046,
    919.0274, 915.0533, 910.9783, 913.6773, 916.6287, 907.9267,
    908.0421, 908.7398, 911.8401, 914.5696, 912.0115, 919.4418,
    917.0436, 920.5495, 917.6138, 907.5037, 908.5145, 919.5846,
    917.6047, 926.8447, 910.6347, 912.8305, 907.7085, 911.6889,
]

for n in xrange(1, 6):
    print 'Looking for %d maxima:' % n
    indexes = find_n_maxima(data, n)
    print indexes
    print ', '.join(str(data[i]) for i in indexes)
    print

Вывод:

Looking for 1 maxima:
[78]
2881.42

Looking for 2 maxima:
[78, 218]
2881.42, 1927.083

Looking for 3 maxima:
[78, 108, 218]
2881.42, 2147.193, 1927.083

Looking for 4 maxima:
[78, 108, 218, 274]
2881.42, 2147.193, 1927.083, 936.7066

Looking for 5 maxima:
[78, 108, 218, 269, 274]
2881.42, 2147.193, 1927.083, 939.6088, 936.7066

Answer 1

Не могли бы вы двигаться по графику, регулярно вычисляя производную, и если она переключается с положительной на отрицательную, вы нашли пик?

Answer 2

Локальными максимумами будет любая точка x, которая имеет более высокое значение y, чем любой из ее левого и правого соседей. Чтобы устранить шум, вы можете установить некоторый порог допуска (например, точка x должна иметь более высокое значение y, чем n ее соседей).

Чтобы не сканировать каждую точку, вы можете использовать тот же подход, но использовать 5 или 10 баллов за раз, чтобы получить приблизительное представление о том, где лежит максимум. Затем вернитесь в эти области для более детального сканирования.

Answer 3

прямой подход, что-то вроде этого:

#include <stdio.h>
#include <stdlib.h>

#define MAXN 100

double smooth(double arr[], int n, int i)
{
        double l,r,smoo;
        l = (i - 1 < 0)?arr[0]:arr[i-1];
        r = (i + 1 >= n)?arr[n-1]:arr[i+1];
        smoo = (l + 2*arr[i] + r)/4;
        return smoo;
}

void findmax(double arr[], int n)
{
        double prev = arr[0];
        int i;
        int goup = 1;

        for(i = 0; i < n; i++)
        {
                double cur = smooth(arr,n,i);
                if(goup) {
                        if(prev > cur && i > 0) {
                                printf("max at %d = %lf\n", i-1, arr[i-1]);
                                goup = 0;
                        }
                } else {
                        if(prev < cur)
                                goup = 1;
                }
                prev = cur;
        }
}

int main()
{
        double arr[MAXN] = {0,0,1,2,3,4,4,3,2,2,3,4,6,8,7,8,6,3,1,0};
        int n = 20, i;

        for(i = 0; i < n; i++)
                printf("%.1lf ",arr[i]);
        printf("\n");

        findmax(arr,n);
        return 0;
}

вывод:

0.0 0.0 1.0 2.0 3.0 4.0 4.0 3.0 2.0 2.0 3.0 4.0 6.0 8.0 7.0 8.0 6.0 3.0 1.0 0.0
max at 6 = 4.000000
max at 14 = 7.000000

1) set state = goup: движение вверх по кривой;
2), если ранее значениебольше текущего, чем был максимум:
распечатайте его и установите для состояния godown
3) в то время как в состоянии Godown дождитесь, пока предыдущий не станет меньше текущего и переключитесь на (1)

, чтобы уменьшить шумиспользовать некоторую функцию сглаживания smooth ()

Answer 4

Я практикую, я нашел то, что хорошо работает, это использовать операцию морфологии расширения для создания расширенной версии вашей функции выборки (точки данных), а затем для определения локального максимума сравнивать расширенную версию с оригинальной и где угодно расширенная версия, равная исходной версии, должна иметь локальный максимум. Это работает очень хорошо, я нахожу с 2D + данными (то есть изображениями), но так как у вас есть 1D данные, может быть проще просто использовать различия между последовательными точками как приближение к производной.

Обратите внимание: если вы используете метод расширения, то структурирующий элемент (размер и форма), который вы используете при расширении, в значительной степени определяет типы пиков, которые вы ищете.

Также, если у вас есть шум в ваших данных, сгладьте их с помощью фильтра нижних частот, например, 1-мерного гауссова, перед поиском.

Более подробную информацию о расширении можно найти здесь:

вот идея, реализованная в matlab: http://www.mathworks.com/matlabcentral/fileexchange/14498-local-maxima-minima

если вы не знаете, что такое дилатация: http://en.wikipedia.org/wiki/Dilation_%28morphology%29

(это очень просто, как только вы это поймете, это очень простое объяснение) http://homepages.inf.ed.ac.uk/rbf/HIPR2/dilate.htm

Answer 5

Вы можете попытаться подогнать сплайн к данным, а затем найти экстремумы сплайна. Поскольку сплайн является кусочно-полиномиальным, точные положения экстремумов можно найти с помощью относительно простых формул.

Answer 6

Другой метод состоит в том, чтобы создать то, что я называю средним числом пешеходных склонов. Я не знаю, есть ли имя для него, но оно выглядит так, ваш набор данных, например, 1000 чисел, вы берете x [n] + x [n..1234567] чисел скажем, 7 чисел, в среднем, первые 3 и последние 3 используйте их, чтобы узнать, будет ли линия над ними идти вверх или вниз.

Когда он падает, вы проходите номер горного пика, после того, как один такой образец продолжает ждать, пока он снова не поднимется. Так что только после того, как вверх, вы отслеживаете первый момент спуска. И повтори это.

Он обнаружит вершины, и в зависимости от длины линии наклона (7) .. 15 .. 33 и т. Д., Вы также удалите шум.

Answer 7

Поскольку другие упоминали производные или сравнение с местными соседями, обычно работает для меня. Если вас беспокоит шум, я могу порекомендовать медианная фильтрация в качестве довольно быстрой и надежной схемы фильтрации. Я использую его и все время обращаю медианную фильтрацию, чтобы подавить шум в акустических датчиках, прекрасно работает.

Answer 8

Вы можете попробовать полосовой фильтр, чтобы подавить шум и упростить надежный выбор этих максимумов.

Смысл полосового (а не низкочастотного) тяготения - почтипостоянная до нуля.Таким образом, самые высокие значения, которые вы найдете в отфильтрованном результате, вероятно, будут самыми ясными пиками.

Конечно, если вы можете определить узкий частотный диапазон для вашего сигнала и применить очень избирательный фильтр, это должно сделать довольнопрактичный простой алгоритм поиска максимума - например, простое сканирование "выборки выше, чем у соседей".

Сложный фильтр может не понадобиться - вам может понадобиться мексиканский вейвлет шляпы в единственном удачно выбранном масштабе.Одна шкала, вероятно, означает, что это больше не вейвлет - просто полосовой FIR-фильтр.

EDIT

Есть асимметричный вейвлет (я забыл название), который, если мексиканская шляпа аналогична косинусу, берет на себя роль синуса.Я упоминаю об этом, поскольку он объединяет полосовую фильтрацию с некой производной - пересечения нуля в результате являются стационарными точками полосовой фильтрованной версии исходного сигнала.сканирование может затем идентифицировать надежные максимумы путем поиска точек пересечения в этом «производном» сигнале.

Answer 9

Знаете ли вы производные данных?Если да, и вы можете символически решить систему, то вы можете найти точки, где производная равна нулю.

Если у вас нет формулы (что, как кажется, имеет место из вашего ОП), тогда выВы можете попробовать отфильтровать некоторый шум, а затем сделать следующее:

Если вы не можете решить символически, тогда вы можете использовать что-то вроде метода Ньютона-Рафсона, чтобы добраться до локальных максимумов и выбрать начальные точки случайным образом издиапазон, чтобы попытаться захватить все максимумы.

Если у вас нет производных данных, вы можете попробовать алгоритм восхождения на гору, который не требует производной, и запустить его в нескольких случайно выбранных точках.Затем вы можете отслеживать точки, которые вы найдете, когда завершается итеративная часть алгоритма восхождения на гору.Это только вероятностно найдет все локальные максимумы, но может оказаться достаточно для ваших целей.

РЕДАКТИРОВАТЬ: учитывая, что 3 пика будут примерно в тех же местах, вы должны попытаться гарантировать, что отправная точка для этих алгоритмовнаходится рядом с этими точками хотя бы в некоторых случаях, когда вы запускаете итерационный алгоритм.